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Abstract 

The  primary  objective  of  this  research  is  to  extend  fully  nonlinear  potential  flow  theory  to 
the  problem  of  breaking  wave  impact  upon  a  vertical  wall.  The  validity  of  this  theory  as  a 
basis  for  fhe  simulation  of  a  deep-water  plunging  breaking  wave^ha^  recently  'beeir 
demonstrated  by  comparison  of  numerical  results  with  experiments  by  Chan.  The 
numerical  approach  used  was  a  mixed  Eulerian-Lagrangian  method  based  upon  application  ^  ..  ^  , 

of  Cauchy’s  integral  theorem  to  the  fluid  domain.,*  This  same  approach  is^appliedrto  the 
wave  impact  problem.  &  the  present  work.— The  experimental  results  of  Chah-on  the  -- 

kinematics  and  dynamics  of  wave  impact  on  a  vertical  wall  provide  a  basis  for  comparison. 

As  a  preliminary  study,  the  simulation  of  a  plunging  breaker  in  a  long  tank  is  repeated  to 
examine  the  difficulties  involved  more  closely  and  to  investigate  the  influence  of  point 
regridding  and  smoothing  routines  upon  the  results.  Next, 'the  a  shorter  numerical  tank  is 
used  with  the  same  wavemaker  input  in  order  to  force  overturning  to  occur  near  the  end  of 
the  tank.  'Several-Difficulties  with  simulation  are  revealed  which  were  not  present  with  the 
longer  tank,  and  techniques  are  developed  to  address  them.  A  simulation  is  finally  achieved 
in  which  the  overturning  wave  peaches  the  end  of  the  tank  before  re-entry  into  the  free 
surface  occurs.  This  allows  ’the  imposition  of  an  impact  condition  at  the  end  of  the  tank. 

-Various techniques  for  simulating  the  impact  process  are  then  investigated. 

The  preliminary  results  obtained  for  the  impact  simulations  show  qualitatively  plausible 
flows  in  some  cases,  with  the  formation  of  upward  and  downward  moving  jets  along  the 
wall  after  impact.  The  computed  impact  pressures  on  the  wall,  however,  are  initially 
negative  and  are  found  to  be  dependent  on  the  number  of  points  used  to  define  the  impact 
region  on  the  wall.  Further  investigation  is  required  to  assess  the  ability  of  this  potential 
flow  model  to  simulate  wave  impact. 
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Chapter  1 
Introduction 

Two  natural  processes  which  are  little  understood  from  a  theoretical  viewpoint  are 
the  evolution  of  steep  and  overturning  water  waves  and  the  impact  of  such  waves  upon 
rigid  structures.  The  development  of  numerical  models  to  gain  insight  into  these  processes 
is  of  great  interest,  especially  in  the  case  of  wave  impact  on  structures. 

The  ability  of  a  fully  nonlinear  potential  flow  model  to  describe  the  motion  of  steep 
and  overturning  gravity  waves  was  confirmed  in  the  work  of  Dommermuth,  Yue,  Lin, 
Rapp,  Chan  and  Melville  [3], in  which  numerical  simulations  based  upon  nonlinear  potential 
theory  reproduced  experimental  measurements  of  surface  elevations  and  fluid  velocities  at 
various  depths  for  a  plunging  breaker  generated  in  a  wave  tank  by  a  piston-type 
wavemaker.  The  numerical  approach  used  was  a  mixed  Eulerian-Lagrangian  scheme 
similar  to  that  developed  by  Vinje  and  Brevig  [9].  The  success  of  a  model  which  ignores 
such  physical  phenomena  as  viscosity  and  surface  tension  also  implies  that  these 
mechanisms  are  relatively  unimportant  for  plunging  breaking  waves  [3],  at  least  far  from 
surface-piercing  boundaries  and  prior  to  re-entry  of  the  overturning  crest  into  the  free 
surface. 

The  kinematics  and  dynamics  of  deep-water  plunging  waves  impacting  on  a  vertical 
wall  were  investigated  experimentally  by  Chan  [1]  using  the  same  wavemaker  and  tank 
setup.  In  the  experiment,  a  computer-controlled  wave  generator  was  first  used  to 
consistently  reproduce  a  plunging  breaker  at  location  in  the  tank  far  from  either  end.  A 
movable  surface  piercing  plate  was  then  positioned  at  various  locations  in  the  observed 
zone  of  "open-water"  wave  breaking  so  that  a  range  of  impact  conditions  on  the  wall  was 
achieved.  Systematic  measurements  of  wall  pressures,  fluid  velocities  and  surface 
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elevations  were  made  for  each  condition  which  allowed  the  correlation  of  the  kinematic  and 
dynamic  properties  of  the  impact  process.  One  of  the  main  observations  made  was  that 
wave  impact  occurs  as  the  wave  front  is  "focused"  toward  a  zone  on  the  wall  termed  the 
"impact  zone",  and  that  air  is  trapped  in  the  impact  process  to  varying  degrees.  It  was 
found  that  the  impulsive  pressures  measured  at  impact  were  highest  when  a  relatively  large 
pocket  of  air  was  trapped  at  impact.  Pressure  oscillations  were  also  measured  at  impact  and 
were  attributed  to  the  dynamics  of  the  trapped  air,  and  it  was  concluded  that  the  impact 
dynamics  could  be  decomposed  into  a  non-oscillating  hydrodynamic  component  and  an 
oscillating  pressure  component  due  to  trapped  air.  The  impact  pressure  characteristics  were 
found  to  be  strongly  dependent  on  the  wall  location  relative  to  the  wave  breaking  location, 
although  significant  variation  in  impact  pressures  was  observed  at  a  given  wall  location  due 
primarily  to  the  random  nature  of  the  trapped  air  dynamics  [2]. 

The  main  purpose  of  the  present  work  is  to  investigate  the  application  of  fully 
nonlinear  potential  flow  theory  to  the  problem  of  a  plunging  breaker  impacting  a  vertical 
wall,  using  the  experimental  results  of  Chan  [1,2]  as  a  basis  for  comparison.  The  starting 
point  for  the  investigation  is  the  same  mathematical  and  numerical  approach  and  computer 
program  used  for  the  plunging  breaker  simulation  in  a  long  tank  [3].  From  this  point  the 
primary  objectives  of  the  research  are  to: 

1.  Extend  the  program  developed  by  Dommermuth  et  al.  [3]  to  permit  simulation 
of  the  impact  of  an  overturning  wave  upon  a  vertical  wall;  and ' 

2.  Investigate  criteria  based  upon  local  flow  characteristics  which  may  be  used 
to  develop  a  robust  rule  to  automatically  control  the  distribution  of 
Lagrangian  points  on  the  free  surface  during  the  simulation. 

The  second  objective  is  important  if  a  general  algorithm  is  to  be  found  which  will 
reduce  or  eliminate  the  need  for  trial-and-error  approaches  to  the  use  of  such  numerical 
techniques  as  regridding  and  smoothing  in  order  to  complete  a  simulation. 
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Chapter  2 

Mathematical  Formulation 

The  assumptions  which  allow  potential  theory  to  be  applied  are  that  the  fluid  flow  is 
irrotational  and  inviscid.  The  fluid  is  also  assumed  to  be  incompressible  and  homogeneous 
with  gravity  as  the  only  body  force  acting.  Surface  tension  is  also  neglected.  The  fluid  is 
confined  to  a  two-dimensional  rectangular  tank  of  length  L{t)  with  a  vertical  piston 
wavemaker  at  one  end  which  undergoes  a  prescribed  time-dependent  horizontal  motion 
with  velocity  (/(/). 

There  are  two  problems  to  examine:  the  case  of  flow  prior  to  impact  with  the  wall 
and  the  case  of  flow  after  impact.  The  mathematical  formulation  for  the  two  cases  differs 
only  in  the  boundary  conditions  used,  and  is  similar  to  the  approach  used  by  Vinje  and 
Brevig  [9]  in  that  it  is  based  upon  Cauchy’s  integral  theorem.  The  equations  below  are 
presented  Li  the  form  given  by  Dommermuth  et  al.  [3].  Note  that  mass,  length  and  time 
units  are  chosen  such  that  the  non-dimensionalized  gravitational  acceleration,  density  and 
tank  depth  have  a  value  of  unity.  This  applies  to  all  equations  and  results  presented 
hereafter. 

A  sketch  of  the  fluid  domain  is  given  in  Figure  2-1  for  the  case  where  wave  impact 
with  the  right  boundary  has  not  yet  occurred. 

The  complex  potential  is  defined  as  p(z,r)  =  <| +  i\yU,v,r),  where  both  <j>  and  \\l 
are  functions  which  satisfy  Laplace’s  equation  and  z  =  x+iy.  If  the  entire  boundary  is 
denoted  as  C(z,t)  =  Bl\j  F  Br^j  Bb,  then  Cauchy’s  integral  theorem  yields 

2tti  |3(£,r)  =  J  ^£d z  2.1 

j  Q  z  C, 


) 
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Bb 

Figure  2- 1 :  Fluid  domain  for  pre-impact  flow, 
where  C  is  inside  the  contour  of  integration  The  boundary  conditions  for  the  pre-impact 
problem  are  given  below: 


V  =  U(r)  (y  +  1 )  on  BL(x.t) 

2.2 

y  =  0  on  Br  ,  Bb 

2.3 

=  on  F(z.r) 

Dr  & 

2.4a 

^  =  0.5  V/2  -  v  -p  on  F(:,t) 

2.4b 

In  the  pre-impact  case  the  pressure  p  on  the  free  surface  is  set  to  zero.  Note  that  that 
D  /Dr  a  d  /9r  +  V< j>  V  is  the  material  derivative  and  that  *  represents  the  complex 
conjugate.  Note  also  that  \  =  u  +  v*  =  (dp/cfc)  (d(3/dz)*,  where  V  is  the  velocity 

magnitude,  u  is  the  horizontal  component  of  velocity  and  v  is  the  vertical  component.  The 
specification  of  the  problem  for  f3(:,r)  and  F(r.r)  is  complete  when  initial  conditions  are 
given  which  correspond  to  a  fluid  at  rest  at  r  =  0 
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As  in  Dommermuth  et  al.  [3],  the  fact  that  $  is  even  with  respect  to  the  bottom 
boundary  Bb  and  y  is  odd  with  respect  to  Bb  is  used  to  reduce  the  number  of  unknowns  by 
the  method  of  images  If  the  contour  of  integration  is  reflected  about  the  bottom  line  v  =  -1 
and  the  bottom  is  then  removed.  Cauchy's  integral  theorem  then  yields 


2rn  p(C.n  =  j 

J c  L;-b 


2.5 


where  C'  =  B L\j  F  u  B  R.  Numerically,  for  a  tank  with  L  >>  1.  the  number  of 
unknowns  is  substantially  reduced  by  the  number  of  points  defined  originally  along  the 
bottom. 

In  case  of  wall  impact,  the  fluid  domain  is  modified  as  shown  in  Figure  2-2. 


*s 

Figure  2-2:  Fluid  domain  for  post-impact  flow 


Now  C(:.r)  =  F  j  ~  BR  |  uf  2  ^  ^2  ^  boundary  conditions  for  the 

post-impact  problem  are  the  same  except  that  BR  is  divided  into  BRl  and  BR~,,  and  F  is 
divided  into  F j  and  Fy  In  equation  (2  4b ).  a  non-zero  pressure  p  along  F-,  would  represent 
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air  being  trapped  at  impact,  in  which  case  an  additional  relationship  between  p  and  the 
volume  of  the  trapped  air  pocket  is  necessary.  For  the  post-impact  case,  equation  (2.5)  is 
the  same  with  C  modified  appropriately.  In  both  cases,  if  £  is  allowed  to  approach  C  from 
the  interior,  a  Fredholm  integral  equation  of  the  second  kind  for  the  unknowns  y  on  F  and 
on  B  is  obtained  by  taking  the  real  part  of  (2.5)  when  £  is  on  F.  Similarly,  a  Fredholm 
integral  equation  of  the  second  kind  for  the  unknowns  is  obtained  by  taking  the  imaginary 
part  of  (2.5)  when  £  is  on  B  [9]. 

At  a  given  value  of  time  r,  the  exact  solution  of  the  boundary  value  problem  is 
approximated  by  the  solution  of  the  system  of  equations  obtained  by  discretization  of  the 
contour  of  integration  in  equation  (2.5)  to  obtain  the  values  of  y  on  the  free  surface  and  ^ 
on  the  tank  boundaries.  The  complex  velocity  potential  P  is  then  known  on  the  free  suface 
so  that  the  velocity  at  the  free  surface  points  can  then  be  computed  from  d|3/dr  =  u  -  iv. 
Using  these  velocities,  numerical  time  integration  of  the  kinematic  and  dynamic  boundary 
conditions  on  the  free  surface,  equations  (2.4a)  and  (2.4b).  gives  the  position  and  velocity 
potential  on  F{z,t  +  Ar)  to  be  used  in  solving  the  BVP  at  time  t  +  At.  The  wavemaker 
velocity  U(t)  is  always  a  known  input  to  the  problem.  The  solution  of  the  BVP  is 
independent  of  time,  all  points  being  considered  Eulerian  points,  while  the  time  integration 
step  treats  the  free  surface  points  as  Lagrangian  points;  hence  the  description  of  solution  as 
a  mixed  Eulerian-Lagrangian  approach. 


Chapter  3 

Numerical  Implementation 


3.1  Discretization,  Solution  and  Time  Integration 

The  discretization  of  equation  (2.5)  is  achieved  by  dividing  the  modified  contour  of 
integration  C  into  linear  segments  or  panels.  The  complex  potential  P  is  assumed  to  vary 
linearly  between  the  endpoints  of  each  panel,  which  are  referred  to  as  nodes.  At  a  given 
node  2^  which  corresponds  to  £  in  (2.5),  integration  along  the  panelized  contour  with  nodes 
Zj,  which  correspond  to  the  variable  of  integration  z  in  (2.5),  yields  one  equation  relating  the 
unknowns  on  the  boundary.  The  real  or  imaginary  part  of  the  resulting  discretized  equation 
is  used  appropriately  to  maintain  the  foim  of  a  Fredholm  integral  equation  of  the  second 
kind,  as  discussed  in  Chapter  2.  When  this  is  repeated  for  all  zk,  an  N  x  N  linear  system  of 
equations  Ax  =  b  for  the  N  unknowns  is  obtained  where  A  is  the  matrix  of  influence 
coefficients.  The  formulas  for  these  influence  coefficients  and  their  asymptotic 
approximations  for  Zj  far  from  zk  are  given  in  Vinje  and  Brevig  [9]. 

In  this  numerical  approach,  the  treatment  of  the  singularity  of  the  solution  at  the 
intersection  point  between  the  rigid  boundaries  and  the  free  surface  is  the  same  as  that 
developed  by  Lin  [5, 6]:  both  <t>  and  \| f  are  considered  known  values  at  the  contact  points. 
This  approach  is  also  used  in  extending  the  program  used  by  Dommermuth  et  al.  [3]  to  the 
case  for  wave  impact  in  Chapter  6,  where  two  new  contact  points  are  considered  to  be 
created  at  the  right  wall  when  the  plunging  crest  tip  impacts  the  wall. 

The  solution  of  the  boundary  value  problem  at  each  time  step  is  performed  using  the 
standard  LINPACK  routines  SGEFA  for  factorization  and  SGESL  for  back-substitution. 
After  solution,  the  complex  potential  p  is  known  at  all  points  on  the  boundary.  A  three- 
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point  finite  difference  formula  is  then  used  to  compute  the  velocities  of  the  Lagrangian 
points  on  the  free  surface,  since  d$/dz  =  u  -  iv.  The  horizontal  and  vertical  velocities  u  and 
v  then  specify  the  right-hand  sides  of  equations  (2.4),  allowing  time  integration  of  these 
equations  to  obtain  new  known  values  of  <j>  and  z  on  F  at  the  next  time  step. 

The  method  of  time  integration  used  is  a  combination  of  a  fourth-order  Runge-  Kutta 
(RK4)  scheme  to  integrate  values  for  three  initial  time  steps  followed  by  a  fourth-order 
Adams-Moulton-Bashforth  (ABM4)  scheme  for  an  arbitrary  number  of  steps.  This 
technique  was  also  used  by  Longuet-Higgins  and  Cokelet  [7],  and  the  linear  stability 
aspects  of  each  scheme  are  addressed  in  Dommermuth  et  al.  [3]. 

3.2  Regridding  and  Smoothing;  Time  Step  Control 

The  program  developed  by  Dommermuth  et  al.  [3],  which  will  henceforth  be  referred 
to  as  the  "original"  program,  uses  two  techniques  to  avoid  or  supress  ’’sawtooth’’ 
instabilities  in  the  evolution  of  the  free  surface  which  have  been  observed  by  several 
investigators  (e  g.,  [7])  using  the  mixed  Eulerian- Lagrangian  approach.  One  technique  is  a 
five-point  smoothing  formula  identical  to  that  used  by  Longuet-Higgins  and  Cokelet  [7], 
which  is  applied  to  values  of  position  and  velocity  potential  at  the  free  surface  points.  The 
other  is  a  regridding  algorithm  [3, 4)  which  uses  quadratic  interpolation  of  position  and 
velocity  potential  at  the  free  surface  points  to  create  a  new  set  of  equally  distributed  points 
using  the  linear  distance  between  the  original  points  as  a  parameter.  Both  smoothing  and 
regridding  tend  to  suppress  the  effects  of  higher-wavenumber  components;  qualitatively, 
smoothing  treats  any  oscillatory  type  of  variation  as  noise  and  simply  finds  an  average 
curve  through  the  input  values,  while  regridding  tends  to  prevent  the  resolution  of  higher 
wavenumbers  by  arresting  the  concentration  of  Lagrangian  points  in  regions  of  higher 
velocity  gradients  and  wavenumbers.  The  cause  of  "sawtooth”  instablities  is  thought  to  be 
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related  in  part  to  possible  growth  of  errors  in  computed  velocities  due  to  small  errors  in 
computed  velocity  potential  in  such  regions  [3].  In  the  original  program,  the  regridding 
alogorithm  is  applied  to  all  points  on  the  free  surface  between  the  wavemaker  and  far  wall. 
As  a  wave  steepens  and  eventually  overturns,  a  certain  degree  of  spatial  resolution  is 
required  to  define  regions  of  high  velocities  and  curvature,  such  as  the  crest  tip.  For  this 
reason,  the  approach  taken  in  the  simulation  of  a  plunging  breaker  using  the  original 
program  is  to  employ  regridding  until  more  resolution  is  required,  at  which  time  smoothing 
is  invoked  in  order  to  allow  the  free  surface  points  to  "cluster"  in  the  overturning  region.  A 
dynamic  time-step  control  algorithm  is  also  used  which  limits  the  time  step  size  so  that  no 
panel  moves  more  than  10%  of  its  length  in  one  time  step  [3]. 

In  Chapter  4,  the  original  program  is  used  to  do  further  simulations  for  a  plunging 
breaker  using  the  same  wavemaker  input  and  tank  length  (L=20)  as  Dommermuth  et  al.  [3]. 
The  purpose  of  these  numerical  experiments  is  to  investigate  the  dependence  of  results  on 
the  particular  times  when  regridding  is  stopped  and  smoothing  is  started.  In  Chapter  5, 
modification  in  the  use  of  smoothing  at  the  free  surface-wall  intersection  point  is  found 
necessary  to  simulate  a  plunging  breaking  wave  just  prior  to  impact  with  a  vertical  wall. 
The  regridding  algorithm  is  also  modified  to  add  points  in  the  curl  region  to  obtain 
sufficient  resolution.  This  trial-and-error  approach  to  the  use  of  smoothing  and  regridding 
is  the  motivation  for  Chapter  7,  which  investigates  the  possibility  of  developing  a  robust 
method  of  free  surface  point  distribution  based  upon  criteria  related  to  local  flow  properties. 

3.3  Conservation  of  Mass  and  Energy 

The  conservation  of  mass  is  is  checked  by  comparing  the  original  fluid  area  (volume) 
with  a  numerical  computation  of  area  based  on  the  following  formula: 
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A  = 


3.1 


In  this  formula,  C  represents  the  original  closed  contour  F  u  Br^j  Bb  and  a 
trapezoidal  rule  is  used  for  the  free  surface  integration. 

The  conservation  of  energy  is  checked  in  one  form  by  comparing  the  computed 
power  input  from  the  wavemaker  with  the  rate  of  change  of  fluid  energy,  represented  in 
continuous  form  as  follows: 

The  left-hand  side  represents  the  power  input  from  the  wavemaker.  The  first  term  in 
brackets  on  the  right  hand  side  represents  the  kinetic  energy  of  the  fluid  and  the  second 
term  is  the  potential  energy  contribution  from  the  free  surface.  The  last  term  represents  the 
rate  of  change  of  potential  energy  due  the  changing  tank  length,  and  in  practice  is 
subtracted  from  both  sides  to  improve  the  accuracy  of  the  estimate  [3],  A  trapezoidal  mle 
is  again  used  for  the  integrations,  and  the  normal  vector  n  =  nxi  +  n} j  at  each  point  on  the 
free  surface  (which  points  outward  from  the  fluid)  is  approximated  by  the  average  of  the 
normal  vectors  to  the  adjacent  panels,  except  at  the  contact  points  where  the  normal  to  the 
first  or  last  free  surface  panel  is  used.  A  second  check  of  energy  conservation  is  obtained 
by  comparing  the  total  fluid  energy  computed  directly  at  each  time  step  with  the  total  work 
input  obtained  by  integrating  in  time  the  left-hand  side  of  equation  (3.2)  using  a  three-point 
integration  routine. 

Note  that  both  (3.1 )  and  (3.2)  are  valid  for  both  the  pre-impact  and  post-impact  cases, 
the  only  difference  being  integration  over  two  free  surface  regions  in  the  post-impact  case. 


3.4  Structure  of  Program 


The  original  program  is  written  in  two  parts:  a  start  program  for  time  equal  to  zero  to 
some  later  value  tn  and  a  restart  program  from  time  tn_{  to  some  time  /y.  The  basic  logic  of 
the  programs  is  given  below: 

1.  INPUT:  known  values  of  x,  v  and  <j>  on  F  and  y  on  B  at  time  step  N. 

2.  RK4  integration  for  steps  N+l  to  N+3.  For  each  sub-step  k=l  to  k=4  in  the 
RK  computation  at  a  given  time  step  do  the  following: 

a.  Compute  value  of  t,  y  and  <J>  on  F  at  sub-step  k  using  the  values  from 
step  (g)  below;  or.  if  k=l,  just  save  values  from  previous  time  step. 

b.  Reposition  nodes  on  wavemaker  and  far  wall  so  that  they  are  evenly 
spaced  based  on  elevation  of  contact  points. 

c.  Compute  the  value  of  wavemaker  velocity  using  linear  interpolation  of 
given  velocity  time  history,  assign  values  of  y  on  BL. 

d.  Compute  the  influence  coefficients  and  set  up  the  A  matrix  and  b 
vector. 

e.  Solve  the  system  for  <j>  on  B  and  y  on  F. 

f.  Compute  the  velocities  u  and  v  at  the  free  surface  points  by  defining 
complex  potential  (3  and  position  z  at  each  point  and  using  three-point 
difference  along  F  (one-sided  for  contact  points  and  centered  for 
interior  points  i  Compute  value  of  D<|>/Dt  in  equation  (2.4a). 

g.  Compute  At  it.  Ar  v  and  A/  D<j)/D t  at  free  surface  points  for  use  in 
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h.  If  k=l,  save  values  of  computed  in  (g)  for  later  use  in  ABM4 
integration.  If  current  time  step  is  N+l,  values  are  saved 
corresponding  to  N,  etc. 

i.  If  k=l,  write  output  for  previous  time  step:  x,  y,  <J>,  and  y  at  all  nodes. 
For  example,  if  current  time  step  is  N+2,  then  values  are  saved  for 
time  step  N+l. 

j.  If  k=l  and  N  is  the  initial  time  step,  call  time  step  control  alogorithm 
to  set  value  of  A t  using  the  velocities  computed  in  step  (f). 

k.  If  k=4,  complete  RK  integration  for  x,  y  and  (j>  on  F  for  current  time 
step.  If  current  step  is  N+3,  then  proceed  to  ABM4  integration. 

3.  ABM4  integration  for  time  steps  N+4  to  N+M,  where  M  is  arbitrary.  For 
each  sub-step  k=l  to  k=2  in  the  ABM  computation  at  a  given  time  step  do  the 
following: 

a.  Compute  value  of  ,r,  y  and  on  F  at  sub-step  k  ("predictor"  step)  using 
the  values  from  step  (g)  below;  or,  if  k=l,  just  save  values  from 
previous  time  step. 

b.  Reposition  nodes  on  wavemaker  and  far  wall  so  that  they  are  evenly 
spaced  based  on  elevation  of  contact  points. 

c.  Compute  the  value  of  wavemaker  velocity  using  linear  interpolation  of 
given  velocity  time  history,  assign  values  of  y  on  BL. 

d.  Compute  the  influence  coefficients  and  set  up  the  A  matrix  and  b 
vector. 

e.  Solve  the  system  for  <(>  on  B  and  y  on  F. 
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f.  Compute  the  velocities  u  and  v  at  the  free  surface  points  by  defining 
complex  potential  (3  and  position  z  at  each  point  and  using  three-point 
difference  along  F  (one-sided  for  contact  points  and  centered  for 
interior  points).  Compute  value  of  D<j)/Dt  in  equation  (2.4a). 

g.  Compute  At  ■  u,At  ■  v  and  At  ■  D<p/D t  for  use  in  (a). 

h.  If  k=l,  write  output  for  previous  time  step:  x,  y,  <f>,  and  y  at  all  nodes. 

For  example,  if  current  time  step  is  N+4,  then  values  are  saved  for 
time  step  N+3.  Note  that  if  the  program  stops  at  step  (4)  below,  then 
the  last  values  saved  are  for  time  step  N+M-l. 

i.  If  k=2,  complete  RK  integration  for  current  time  step  ("corrector" 
step).  If  current  time  step  is  N+M,  then  proceed  to  step  (4). 

4.  If  maximum  number  of  time  steps  or  maximum  value  of  simulation  time  is 
exceeded,  then  STOP. 

5.  Regrid  or  smooth  values  of  x,  y  and  <j>  for  free  surface  points. 

6.  Call  time  step  control  algorithm  to  set  new  value  of  At  using  the  most  recent 
computed  velocities  of  free  surface  points  in  step  (2)(f),  k=2. 

7.  Return  to  start  of  RK  integration  loop  at  step  (2)  with  N=N+M. 

Note  that  the  integer  M  determines  the  frequency  at  which  the  RK  integration  is 
restarted.  In  principle,  it  is  necessary  to  return  to  the  RK  integration  routine  to  provide  new 
initial  values  for  the  ABM  routine  if  either  (1)  the  size  of  the  time  step  changes,  or  (2) 
regridding  or  smoothing  free  surface  points,  or  any  other  operation  which  changes  the 
identity  of  the  points  being  "tracked"  in  the  integration,  is  used.  In  the  original  programs 
and  all  modified  versions  discussed  in  this  work  the  time  step  control  routine  as  well  as  any 
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smoothing,  regridding  or  other  routines  which  change  the  identity  of  free  surface  points  are 
called  every  M*11  time  step  even  though  the  specific  routines  and  the  sequence  of  calling 
them  may  differ  from  one  version  to  another.  The  RK  routine  is  therefore  also  restarted  at 
that  time  step. 

Using  the  original  start  program,  regridding  and  smoothing  of  the  free  surface  points 
occurs  every  time  step  just  before  the  final  solution  and  output  is  completed  for  that 
time  step.  As  noted  above,  if  the  program  stops  at  step  (4)  at  current  time  step  N+M,  the 
last  values  saved  are  for  time  step  N+M-l.  The  values  of  <|>,  x  and  y  on  the  free  surface  at 
time  step  N+M-l  have  not  yet  been  smoothed  or  regridded.  When  the  original  restart 
program  is  then  used,  the  values  at  time  step  N+M- 1  are  read  and  smoothing  or  regridding 
is  immediately  applied  to  the  free  surface  point  values.  The  program  then  proceeds  as  in 
the  start  program.  Because  the  regridding  or  smoothing  occurrence  is  shifted  back  one  time 
step  on  restart,  solution  values  ai  a  given  time  t  using  the  start  program  alone,  for  instance, 
may  differ  slightly  from  values  obtained  using  the  restart  program  one  or  more  times  (M  is 
assumed  the  same).  Also,  there  will  not  necessarily  be  an  exact  matching  of  time  and  time 
step  values  between  runs  which  restart  at  different  points,  since  computed  velocities  may 
also  be  slightly  different  and  because  the  time  step  control  routine  computes  At  based  upon 
velocities  and  positions  of  free  surface  points. 

One  error  in  both  original  programs  is  that  an  initial  value  of  At  is  not  input;  therefore 
at  step  (2)(g),  for  k=l,  the  fust  sub-step  in  the  ABM  integration  incorrectly  uses  A/=0.  In 
step  (2Xj),  the  time  step  control  routine  is  called  and  a  value  of  At  is  computed  and  there  are 
no  further  problems.  This  error  is  corrected  by  switching  steps  (2)(g)  and  (2)(j).  Since  this 
error  only  occurs  at  the  initial  step  of  the  start  or  restart  programs,  the  effect  of  the  error  is 
probably  not  significant;  the  results  using  the  start  program  alone  would  be  essentially 
unaffected  since  at  the  initial  time  step  the  values  computed  for  u,  v  and  Cty/Dr  in  step 
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(2Xg)  are  near  zero  at  all  free  surface  points.  Note  also  in  the  original  programs  that  the 
time  step  control  routine  is  called  in  (6)  after  the  final  integration  step  (3)(i)  for  the  values 
of  $,  x  and  y  on  the  free  surface,  but  before  free  surface  point  velocities  are  computed  based 
upon  these  final  values  after  return  to  the  start  of  the  RK  integration  loop  in  step  (2)(f)  with 
k=l.  Thus  the  positions  of  the  free  surface  points  used  in  the  time  step  control  routine  are 
the  current  positions,  but  the  velocities  used  are  not  current.  Although  this  is  not  an  error,  it 
is  more  consistent  to  call  the  time  step  control  routine  after  the  current  velocities  are 
computed. 

Modified  versions  of  the  original  start  and  restart  programs  were  created  which  allow 
for  a  "consistent"  restart.  Using  the  above  example,  this  means  that  the  last  values  saved 
when  either  program  stops  are  those  for  time  step  N+M,  after  any  regridding  or  smoothing 
has  occurred,  rather  than  for  N+M-l.  On  restart,  these  values  from  step  N+M  are  used  as 
input  and  the  solution  at  time  step  N+M  is  duplicated  so  that  the  program  yields  results  as  if 
no  restart  had  been  necessary.  This  is  desirable  when  program  run  time  is  limited  or  when 
it  is  necessary  to  restart  a  run  at  an  earlier  time  step  but  using  different  input  parameters. 

The  error  mentioned  above  as  well  as  the  inconsistency  in  use  of  the  time  step  control 
routine  are  also  eliminated.  The  logic  of  the  modified  programs  is  given  below: 

1 .  INPUT:  known  values  of  x,  y  and  <J>  on  F  and  on  B  at  time  step  N. 

2.  RK4  integration  for  steps  N+l  to  N+3.  For  each  sub-step  k=l  to  k=4  in  the 
RK  computation  at  a  given  time  step  do  the  following: 

a.  If  k=l  and  time  step  N=N+M  from  step  (4),  regrid  or  smooth  values  of 
x,  y  and  4>  for  free  surface  points. 

b.  Compute  value  of  v.  v  and  <()  on  F  at  sub-step  k  using  the  values  from 
step  (k)  below;  or.  if  k=l,  just  save  values  from  previous  time  step. 
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c.  Reposition  nodes  on  wavemaker  and  far  wall  so  that  they  are  evenly 
spaced  based  on  elevation  of  contact  points. 

d.  Compute  the  value  of  wavemaker  velocity  using  linear  interpolation  of 
given  velocity  time  history,  assign  values  of  vj/  on  BL. 

e.  Compute  the  influence  coefficients  and  set  up  the  A  matrix  and  b 
vector. 

f.  Solve  the  system  for  4»  on  B  and  y  on  F. 

g.  Compute  the  velocities  u  and  v  at  the  free  surface  points  by  defining 
complex  potential  P  and  position  z  at  each  point  and  using  three-point 
difference  along  F  (one-sided  for  contact  points  and  centered  for 
interior  points).  Compute  value  of  Ety/Dt  in  equation  (2.4a). 

h.  If  k=l,  write  output  for  previous  time  step:  v.  y.  and  y  at  all  nodes. 
For  example,  if  current  time  step  is  N+2,  then  values  are  saved  for 
time  step  N+l. 

i.  Call  time  step  control  algorithm  to  set  new  value  of  A t  using  the  most 
recent  computed  velocities  of  free  surface  points  at  step  (g). 

j.  If  maximum  number  of  time  steps  or  maximum  value  of  simulation 
time  is  exceeded,  then  STOP. 

k.  Compute  A t  u.  Ar  v  and  A t  Cty/D t  at  free  surface  points  for  use  in 
(b). 

l.  If  k=l,  save  values  of  computed  in  (k)  for  later  use  in  ABM4 
integration.  If  current  time  step  is  N+l,  values  are  saved 
corresponding  to  N.  etc. 
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m.  If  k=4,  complete  RK  integration  for  x,  y  and  $  on  F  for  current  time 
step.  If  current  step  is  N+3,  then  proceed  to  ABM4  integration. 

3.  ABM4  integration  for  time  steps  N+4  to  N+M,  where  M  is  arbitrary.  For 
each  sub-step  k=i  to  k=2  in  the  ABM  computation  at  a  given  time  step  do  the 
following: 

a.  Compute  value  of  x,  y  and  $  on  F  at  sub-step  k  ("predictor"  step)  using 
the  values  from  step  (g)  below  ;  or,  if  k=l,  just  save  values  from 
previous  time  step. 

b.  Reposition  nodes  on  wavemaker  and  far  wall  so  that  they  are  evenly 
spaced  based  on  elevation  of  contact  points. 

c.  Compute  the  value  of  wavemaker  velocity  using  linear  interpolation  of 
given  velocity  time  history,  assign  values  of  y  on  BL. 

d.  Compute  the  influence  coefficients  and  set  up  the  A  matrix  and  b 
vector. 

e.  Solve  the  system  for  $  on  B  and  y  on  F. 

f.  Compute  the  velocities  u  and  v  at  the  free  surface  points  by  defining 
complex  potential  P  and  position  z  at  each  point  and  using  three-point 
difference  along  F  (one-sided  for  contact  points  and  centered  for 
interior  points).  Compute  value  of  D<J>/Dt  in  equation  (2.4a). 

g.  Compute  A r  u,  A/  v  and  A t  Dty/Dt  for  use  in  (a). 

h.  If  k=  l ,  write  output  for  previous  time  step:  x,  y,  4»,  and  y  at  all  nodes. 
For  example,  if  current  time  step  is  N+4,  then  values  are  saved  for 
time  step  N+3.  Note  that  if  the  program  stops  at  step  (4)  below,  then 
the  last  values  saved  are  for  time  step  N+M- 1 . 


i.  If  k=2,  complete  RK  integration  for  current  time  step  (  corrector 
step).  If  current  time  step  is  N+M,  then  proceed  to  step  (4). 

4.  Return  to  start  of  RK  integration  loop  at  step  (2)  with  N=N+M. 

3.5  Wavemaker  Input 

The  only  input  to  the  original  programs  used  by  Dommermuth  et  al,  [3]  is  the 
velocity  time  history  of  the  vertical  wavemaker.  This  time  history  was  obtained  by  taking  a 
centered  finite  difference  of  the  measured  wavemaker  displacement  in  the  experiment  and 

applying  Fourier  analysis  to  the  data  to  obtain  the  following  Fourier  cosine  series: 

72 

U(r)=^  Uncos((ont-en)  3.3 

rt— l 

A  plot  of  the  velocity  time  history  is  shown  in  Figure  3-1.  The  values  of  Un,  a>n  and 
9n  are  given  in  Appendix  B  of  [3]. 

The  same  input  is  used  in  all  simulations  in  the  present  work,  and  corresponds  to  the 
wavemaker  input  used  by  Chan  in  the  experiments  which  are  used  for  comparison  with 
simulated  results. 

The  electrical  input  signal  to  the  wavemaker  in  the  experiment  is  plotted  in  Figure  1 
of  [2].  It  should  be  noted  that  the  wavemaker  veloctity  input  used  in  the  simulations  shown 
in  Figure  3-1  starts  about  4.3  s  later  than  the  input  wave  signal  to  the  wavemaker  at  a  point 
when  this  signal  first  becomes  non -zero  This  difference  corresponds  to  a  non-dimensional 
time  difference  of  about  17.0.  When  any  references  are  made  in  this  work  to  time  values  in 
the  actual  experiment,  the  values  given  will  reflect  a  correction  of  -17.0  so  that  they  will 
compare  directly  to  the  simulated  time  values. 
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3.6  Convergence  Tests 

The  accuracy  of  the  boundary  value  problem  solution  step  in  the  original  programs  is 
evaluated  in  Dommermuth  et  al.  [3]  using  exact  deep-water  Stokes  waves  for  comparison 
with  the  program  results.  The  convergence  of  relative  error  with  panel  size  is  shown  to  be 
quadratic,  and  for  wave  steepnesses  £  less  than  about  0.2,  a  density  of  40  panels  per 
wavelength  is  found  to  yield  a  relative  error  in  free  surface  velocity  of  approximately  0.5%. 
The  modified  programs  for  the  pre-impact  cases  do  not  alter  the  setup  of  the  system  of 
equations,  so  the  same  convergence  is  expected.  The  only  difference  between  the  original 
and  modified  programs  in  the  solution  step  is  that  the  modified  programs  use  the  standard 
UNPACK  routines  SGEFA  and  SGESL  to  solve  the  system  of  equations.  The  original 
program  also  uses  a  vectorized  Gaussian  elimination  routine  with  partial  pivoting,  and  there 
is  negligible  difference  between  results  using  the  two  solvers. 

The  overall  accuracy  of  the  time  simulation  is  also  evaluated  for  the  original 
programs  in  Dommermuth  et  al.  [3],  The  convergence  of  free  surface  elevation  with 
decreasing  panel  size  and  time  step  size  is  demonstrated  by  comparing  computed  free 
surface  elevations  with  and  without  regridding  applied  to  the  free  surface  points  with  the 
elevation  measured  in  the  experiment  at  a  time  midway  to  the  wave  breaking  event.  To 
check  the  validity  the  program  modifications  described  above,  the  same  convergence  test  is 
performed  using  the  modified  programs,  and  in  addition  two  additional  cases  of  smoothing 
only  and  regridding  plus  smoothing  are  considered.  For  al!  runs  M  =  15,  where  M  is 
defined  in  Section  3.4.  The  tank  length  used  for  the  convergence  test  is  L  =  8,  compared  to 
L  =  20  used  in  the  full  simulation  for  the  plunging  breaker.  The  position  along  the  tank 
length  is  x  =  3.17,  and  the  time  from  the  start  of  the  wavemaker  is  t  =  25. The  results  are 
presented  in  Table  3.1  below  for  an  initially  even  distribution  of  panels  along  the 
wavemaker,  free  surface  and  far  wall.  N  represents  the  total  number  of  panels. 


TABLE  3.1:  EVEN  PANEL  DISTRIBUTION 

N 

At 

No  Regrid 

Regrid 

Smooth 

_ 

Reg.+Sm. 

100 

.100 

-.06033 

-.06011 

150 

.075 

-.06519 

-.06547 

-.06496 

200 

.075 

-.06520 

-.06523 

-.06524 

250 

.050 

-.06574 

-.06557 

-.06576 

-.06561 

For  comparison,  another  test  is  run  where  the  number  of  panels  on  the  wavemaker 
and  wall  is  fixed  at  25  and  the  total  number  of  panels  in  each  case  is  the  same  as  in  Table  1. 
In  other  words,  the  width  of  panels  on  the  free  surface  is  larger  and  the  number  of  free 
surface  panels  smaller  than  in  the  first  test  for  a  given  total  number  of  panels,  except  in  the 
last  case  of  250  panels  which  again  are  evenly  distributed.  The  results  are  presented  in 
Table  3.2  below. 


TABLE  3.2:  UNEVEN  PANEL  DISTRIBUTION 

N 

At 

No  Regrid 

Regrid 

Smooth 

Reg.+Sm. 

100 

.100 

-.05507 

-.04244 

-.05075 

-.04148 

150 

.075 

-.06433 

-.06365 

-.06414 

-.06333 

200 

.075 

-.06512 

-.06528 

-.06521 

-.06526 

250 

.050 

-.06574 

-.06557 

-.06576 

-06561 

The  measured  value  in  the  experiment  was  -.0670. 


In  both  Table  3.1  and  Table  3.2  convergence  is  demonstrated  with  increasing  total 
number  of  panels  and  decreasing  time  step  size.  For  the  case  of  250  panels,  which  has  even 
initial  panel  distribution,  the  values  have  converged  to  within  1 .0%  for  all  combinations  of 
regridding  and  smoothing  based  on  comparison  with  the  case  of  200  panels  in  each  table. 
Note  also  that  for  the  case  of  250  panels  the  run  which  uses  no  smoothing  or  regridding  and 
that  which  uses  smoothing  alone  give  results  closest  to  the  experimental  value.  This 
suggests  that  the  use  of  regridding  may  cause  a  greater  loss  of  fidelity  in  the  numerical 
solution.  In  the  cases  having  uneven  initial  panel  distribution,  the  computed  values  for  100 
and  150  total  panels  are  farther  from  the  experimental  value  than  those  for  the 
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corresponding  cases  with  even  panel  distribution.  In  the  case  of  250  panels,  the  number  of 
panels  per  wave  is  about  40  based  on  a  linear  wave  having  a  frequency  of  to  =  2,  the 
frequency  below  which  most  of  the  input  energy  is  concentrated  based  on  the  harmonic 
content  of  the  wavemaker  velocity  input  spectrum  [3].  This  is  roughly  the  panel  density 
used  in  all  the  simulations  described  in  later  chapters  for  both  the  long  and  short  tank  cases. 

In  the  case  of  even  panel  distribution,  the  first  two  columns  of  Table  3.1  should  be 
identical  to  Table  2  in  [3];  however,  the  convergence  shown  in  Table  3.1  is  not  quite  as  fast. 
As  a  check,  the  original  program  used  in  [3]  was  also  run  using  regridding  alone,  and  the 
results  for  free  surface  elevation,  velocity  potential  and  stream  function  agreed  with  the 
modified  program  up  to  at  least  five  significant  figures.  One  possible  explanation  for  the 
discrepancy  could  then  be  that  a  different  method  of  interpolation  of  the  program  output  in 
time  and  in  x  may  have  been  used  in  [3]  to  obtain  values  at  the  desired  point.  The  method 
used  here  is  a  linear  interpolation  first  in  time  followed  by  a  linear  interpolation  in  x. 

3.6.1  Typical  Computing  Times 

All  simulations  discussed  in  this  work  were  performed  on  a  Cray  2  supercomputer. 
As  mentioned  in  Dommermuth  et  al.  [3],  about  80%  of  the  computational  time  required  for 
the  solution  of  the  system  of  equations  for  the  boundary  value  problem  is  devoted  to  the 
assembly  of  the  matrix  of  influence  coefficients  and  about  20%  is  devoted  to  the  LU 
decomposition  of  the  matrix  solution.  In  Chapter  4,  where  the  simulations  involve  549 

t 

unknowns  for  tank  length  L  =  20,  a  typical  run  (e.g.,  run  A3C)requires  about  21  hours  for 
4200  time  steps.  In  Chapter  5,  where  the  number  of  unknowns  is  reduced  to  349  for  the 
shorter  tank  lengths,  a  typical  run  requires  about  8  hours  for  3600  time  steps  (e.g.,  run  C4). 
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Chapter  4 

Simulation  of  a  Plunging  Breaker  in  a  Long  Tank 

The  simulation  of  a  plunging  breaker  performed  by  Dommermuth  et  al.  [3]  involved 
use  of  regridding  every  15  time  steps  (corresponding  to  M  =  15  in  the  discussion  of 
program  structure  in  Chapter  3)  for  the  initial  phase  (about  3000  time  steps)  and  smoothing 
every  5  steps  for  the  last  phase  (about  1000  steps).  Determination  of  the  time  at  which 
regridding  stopped  and  smoothing  started  was  somewhat  difficult  and  a  matter  of  trial  and 
error;  in  fact,  an  intermediate  phase  during  which  both  regridding  and  smoothing  were 
applied  was  required  in  order  to  allow  the  simulation  to  progress  to  overturning. 

The  purpose  of  this  chapter  is  to  repeat  the  simulation  of  the  plunging  breaker  in  a 
tank  of  L  =  20  (non-dimensionalized  on  a  depth  h  of  0.6  m)  using  the  modified  programs  to 
illustrate  the  difficulties  encountered  for  this  case,  and  also  to  investigate  the  dependence  of 
results  on  various  combinations  of  regridding,  regridding  plus  smoothing  and  smoothing 
alone.  Note  that  the  total  number  of  panels  on  the  free  surface  is  500  and  that  there  are  25 
panels  on  both  the  wavemaker  and  right  wall.  For  L  =  20  and  h  =  1  this  gives  an  initial 
panel  size  of  0.04. 

In  this  and  following  chapters,  the  point  of  run  failure  refers  to  the  point  at  which  the 
computation  breaks  down.  This  usually  coincides  with  the  simulation  becoming  physically 
invalid,  such  as  when  an  overturning  crest  re-enters  the  free  surface  or  if  the  free  surface 
intersects  itself. 

Note  that  plots  of  free  surface  profiles  for  the  runs  referred  to  below  are  located  in 
Appendix  A.  The  run  names  start  with  a  letter  which  corresponds  to  the  appropriate 
appendix.  This  same  convention  is  used  for  all  runs  in  later  chapters  as  well.  Each  free 
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surface  profile  plot  in  the  appendix  has  a  heading  at  the  top  of  the  page  which  gives  the  run 
name,  time  step  number,  and  time  value  at  that  time  step.  Numbers  in  parentheses 
immediately  following  time  values  in  the  text  are  the  corresponding  number  of  time  steps  in 
the  simulation. 

4.1  Results  Using  Modified  Code 

4.1.1  Summary  of  Runs  and  Failure  Types 

Run  A1  is  made  with  regridding  alone  applied  every  15  steps.  Run  A1  fails  at  time  t 
=  44.88  (2149).  The  plot  of  free  surface  profile  from  x  =  0  -  20  at  step  2049,  100  steps 
before  failure,  shows  two  predominant  peaks;  the  first  peak,  nearest  the  right  wall  around  x 
=  9.5  is  the  steepest  (note  the  horizontal  and  vertical  scale)  and  is  followed  by  a  second 
peak  around  x  =  6.6.  The  failure  occurs  at  the  first  peak  as  the  two  panels  at  the  top  of  the 
peak  begin  to  overturn  but  then  collapse  into  one  another  as  shown  in  the  detailed  plots.  It 
is  the  second  peak  noted  above  which,  in  successful  simulations  such  as  runs  A3A  -  A3D 
below,  eventually  steepens  and  overturns  after  t  =  51 . 

Run  A2  is  a  restart  of  run  A1  at  r  =  42. 15  (1770)  with  smoothing  alone  applied  every 
5  steps.  Failure  occurs  at  t  =  43.54  (3575).  Again,  the  failure  is  at  the  first  peak  around  x  = 
9.1.  The  nature  of  the  failure  is  different,  though,  from  run  A1  because  the  free  surface 
points  have  become  more  dense  at  the  top  of  the  first  peak  since  regridding  is  no  longer 
applied.  The  crest  overturns,  forming  a  thin  cusp  which  eventually  self- intersects.  This 
failure  suggests  that  in  the  actual  experiment  the  first  crest  may  itself  have  been  a  weakly 
overturning  or  "spilling"  breaking  wave  which  may  not  be  accurately  simulated  by  the 
present  method.  The  failure  in  run  A1  appears  to  be  a  manifestation  of  a  breaking  tendency 
in  the  absence  of  adequate  resolution  of  the  free  surface.  The  a  priori  knowledge  that  the 
second  peak  is  the  one  of  interest  in  this  case  necessitates  some  means  of  suppressing  the 
tendency  of  the  first  peak  to  break  long  enough  for  the  second  to  develop  to  overturning. 
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Run  A3  restarts  run  A1  at  t  -  42.01  (1755)  with  regridding  plus  smoothing  applied 
every  5  steps.  This  run  is  arbitrarily  stopped  at  t  =  54.29  (3260)  since  the  time  when 
overturning  of  the  second  peak  is  expected  to  occur  has  been  exceeded.  In  this  case 
overturning  of  both  peaks  is  suppressed  by  the  combination  of  regridding  and  smoothing. 

Runs  A3A,  A3B,  A3C  and  A3D  are  restarts  of  run  A3  at  the  following  respective 
times  with  smoothing  every  5  steps:  45.98  (2155),  47.04  (2225),  48.01  (2285)  and  49.02 
(2360).  These  runs  successfully  continue  to  overturning  and  plunging  of  the  second  peak 
and  are  summarized  in  the  next  section.  A  sequence  of  free  surface  profiles  showing  the 
overturning  process  is  plotted  for  each  run  in  Appendix  A.3. 

Runs  A3E,  A3F  and  A3G  are  restarts  of  run  A3  at  the  following  respective  times 
using  smoothing  only  every  5  steps:  44.03  (1970),  44.60  (2030)  and  45.00  (2070).  In  these 
cases,  the  breaking  tendency  of  the  first  peak  has  not  been  successfully  suppressed  and 
failure  occurs  in  a  way  similar  to  run  A2.  The  breaking  of  the  first  peak  occurs  later  than  in 
run  A2  and  the  peak  has  a  smaller  amplitude.  In  these  three  cases,  the  later  smoothing 
alone  is  started,  the  later  breaking  occurs. 

Run  A4  restarts  run  A1  at  an  earlier  time  than  run  A3,  at  t  =  40.09  (1575). 
Regridding  plus  smoothing  is  again  applied  every  5  steps.  This  run  is  arbitrarily  stopped  at 
r  =  50.93  (2575). 

Runs  A4A,  A4B,  A4C  and  A4D  are  restarts  of  run  A4  at  the  following  respective 
times  with  smoothing  every  5  steps:  45.99  (2140),  46.99  (2205),  48.03  (2270)  and  49.03 
(2345).  These  runs  successfully  continue  to  overturning  and  plunging  of  the  second  peak 
and  are  summarized  in  the  next  section.  Because  of  their  similarity  to  runs  A3,  A3A,  A3B, 
A3C  and  A3D  ,  runs  A4  through  A4D  are  not  plotted. 


-32- 

4.1.2  Comparison  of  Successful  Runs 

The  following  table  summarizes  the  time  and  position  of  crest  tip  contact  with  the 
free  surface  (re-entry  point)  for  runs  A3A  -  A3D  and  A4A  -  A4D. 


Trs  is  the  time  when  regridding  plus  smoothing  is  started,  Ts  is  the  time  when  smoothing 
alone  is  started,  Tc  is  the  contact  time,  Xc  and  Yc  are  the  contact  position,  At  is  the  time 
step  size  at  contact  and  STEP  is  the  time  step  number.  The  contact  points  above  are 
determined  by  finding  the  first  time  step  at  which  a  point  on  the  crest  tip  passes  below  a  line 
segment  connecting  two  free  surface  points. 

There  is  little  difference  in  the  contact  time  and  position  between  corresponding  runs 
(i.e.,  A3  A  and  A4A).  This  indicates  that  the  time  at  which  regridding  plus  smoothing  is 
invoked,  Trs,  has  the  least  effect  on  the  results.  The  greatest  dependence  is  on  the  time  at 
which  smoothing  only  is  started,  Ts.  The  latest  re-entry  time  at  the  greatest  distance  from 

the  wavemaker  occurs  when  Ts  is  a  minimum  and  there  is  a  consistent  trend  toward  re- 

« 

entry  occurring  earlier  and  closer  to  the  wavemaker  as  Ts  increases,  although  the  difference 
in  re-entry  time  and  position  between  runs  A3B  and  A3C  or  between  A4B  and  A4C  is 
relatively  small. 


In  the  actual  experunent,  contact  of  the  crest  tip  with  the  free  surface  occurred  at 
approximately  Tc  =  52.2  <12.9  s)  and  Xc  =  12.1  (7.25  m),  and  in  the  simulation  of 
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Dommermuth  et  al.  [3]  contact  occurred  at  approximately  Tc  —  52  and  Xc  =  12.2.  As  in  the 
original  simulation,  all  of  the  above  simulations  yield  re-entry  of  the  plunging  breaker  at  an 
earlier  time  than  the  experiment. 


Examination  of  quantities  at  free  surface  points  such  as  velocity,  acceleration  and 
energy  shows  that  just  prior  to  contact  the  maximum  computed  velocity  magnitude  shifts 
from  values  on  the  order  of  1 .0  in  the  curl  of  the  wave  to  values  an  order  of  magnitude 
higher  at  the  crest  tip  or  cusp,  suggesting  that  computational  error  begins  to  dominate  at  the 
tip  where  the  panel  size  is  about  10  times  smaller  than  the  original  panel  length.  Computed 
maximum  acceleration  magnitudes  are  on  the  order  of  5  -  10  g  in  the  wave  curl  just  prior  to 
this  discontinuity  in  velocity;  Dommermuth  et  al.  [3]  observed  a  maximum  acceleration  of 
about  6  g  in  the  curl  of  the  wave  in  their  simulation  just  before  re-entry. 

For  runs  A3A  -  A3D  the  table  below  gives  energy-related  quantities  at  the  time  step 
just  before  the  discontinuity  in  computed  velocity  occurs. 


TABLE  3.2 

RUN 

STEP 

Vm 

Ec 

-Wc/Ec 

A3A 

4558 

1.18 

.04173 

HESM 

A3B 

4416 

1.04 

.03922 

A3C 

4215 

1.06 

.03974 

1.010 

A3D 

4419 

1.31 

_ 

.04341 

0.924 

Vm  is  the  maximum  velocity  magnitude  at  contact,  Ec  is  the  total  fluid  energy  at  contact 
and  Wc  is  the  total  work  input  to  the  fluid  from  the  wavemaker  at  contact.  The  ratio 
-Wc/Ec  is  a  measure  of  conservation  of  energy  and  will  equal  to  unity  if  the  the  computed 
work  input  and  fluid  energy  have  the  same  magnitude  (Ec  +  Wc  =  0).  In  all  cases  above  the 
computed  work  input  to  the  fluid  Wc  =  .04013,  which  means  that  the  variation  in  total 
energy  between  runs  is  due  only  to  the  variation  in  free  surface  position  and  velocity. 
There  is  a  direct  correspondence  between  the  maximum  velocity  in  the  curl  of  the  wave  at 
contact  and  the  total  energy  at  contact.  Although  the  energy  conservation  check  is  only 
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relative,  since  both  work  and  fluid  energy  quantities  are  derived  from  the  simulation,  runs 
A3B  and  A3C  seem  to  be  the  most  consistent  from  a  global  energy  conservation  standpoint. 
It  should  be  noted  that  for  all  the  above  simulations  the  area  or  "mass"  of  the  fluid  domain 
varies  less  that  0.1%  from  the  initial  area. 

4.2  Summary 

In  repeating  the  simulation  of  the  plunging  breaker  for  a  tank  length  of  L  =  20,  a 
"weak  breaking"  event  appears  to  occur  at  the  peak  leading  the  peak  which  eventually 
becomes  the  plunging  breaker.  Regridding  alone  every  15  steps  fails  in  the  vicinity  of  this 
event  as  the  first  peak  steepens,  although  there  is  not  sufficient  resolution  to  clearly  show  a 
breaking  phenomenon.  When  regridding  is  stopped  and  smoothing  alone  started  at  a  time 
before  the  weak  breaking  event,  the  Lagrangian  points  become  more  dense  as  the  first  peak 
steepens  and  the  breaking  event  is  more  clearly  resolved  into  a  thin  plunging  cusp.  The 
simulation  then  fails  when  the  cusp  self-intersects. 

An  intermediate  application  of  regridding  plus  smoothing  every  5  steps  until  the  time 
of  the  weak  breaking  event  has  passed  followed  by  smoothing  alone  every  5  steps  allows 
the  second  peak  to  develop  to  overturning.  The  time  and  position  at  which  re-entry  occurs 
is  seen  to  be  most  dependent  upon  the  time  at  which  the  regridding  plus  smoothing 
application  is  stopped  and  smoothing  only  begins;  as  this  time  Ts  is  decreased,  the  time  of 
re-entry  is  observed  to  increase.  There  is  a  minimum  value  of  Ts,  however,  below  which 
the  tendency  of  the  first  peak  to  break  is  not  sufficiently  suppressed  and  the  simulation  fails. 

Runs  A3  A,  A3B,  A4A  and  A4B  are  very  close  to  one  another  in  time  and  position  of 
crest  tip  contact  with  the  free  surface;  runs  A3A  and  A3B  also  show  the  best  global 
(relative)  energy  conservation  and  show  a  net  loss  of  fluid  energy  compared  to  work  input 
of  one  to  two  percent.  The  differences  in  final  fluid  energy  at  contact  can  be  attributed 


r 
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primarily  to  variations  in  free  surface  velocities  and  position  which  occur  with  the  changing 
duration  of  the  regridding  plus  smoothing  phase  between  runs. 

When  tank  length  is  shortened  and  the  wavemaker  input  remains  the  same,  similar 


difficulties  in  simulation  are  encountered  which  will  be  described  in  the  next  chapter. 
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Chapter  5 

Simulation  of  a  Plunging  Wave  Near  a  Wall:  Pre-Impact  Flow 

In  this  chapter  the  same  wavemaker  input  is  used  with  a  shorter  tank  length  in  an 
attempt  to  simulate  a  plunging  breaking  wave  near  the  wall  at  the  end  of  the  tank.  The 
initial  choice  of  tank  length.  L  =  1 1 .667.  is  identical  to  one  of  the  three  tank  lengths  used  by 
Chan  in  investigating  the  influence  of  wall  position  on  wave  impact  pressure  [1].  The 
shortest  length  of  the  three  is  chosen  because  of  the  observed  tendency  of  the  simulated 
wave  to  break  sooner  than  in  the  actual  experiment  for  the  plunging  breaker  in  the  longer 
tank,  coupled  with  the  fact  that  in  the  actual  experiment  for  wave  impact  this  tank  length 
corresponded  to  the  least  developed  breaking  prior  to  impact.  Anticipating  earlier  breaking 
than  experiment  for  the  short  tank  cases  as  well,  it  was  thought  that  choice  of  the  shortest 
might  still  permit  simulation  of  a  wave  that  has  started  to  overturn  but  has  not  yet  contacted 
the  free  surface  prior  to  wall  impact.  In  fact,  as  described  below,  re-entry  is  imminent  prior 
to  reaching  the  end  of  the  tank  for  L  =  1 1.667,  and  a  series  of  experiments  with  different 
tank  lengths  is  required  to  obtain  a  simulation  wherein  the  overturning  crest  tip  is 
approximately  horizontal  as  it  approaches  the  end  of  the  tank  so  that  eventual  simulation  of 
impact  with  the  wall  prior  to  re-entry  will  be  possible.  The  case  where  L  =  1 1.60  is  finally 
chosen  to  continue  the  impact  study  in  the  next  chapter.  In  all  cases,  the  number  of  panels 
on  the  free  surface  is  300  and  the  number  of  panels  on  the  wavemaker  and  wall  is  25.  The 
initial  free  surface  panel  size  is  slightly  smaller  than  the  runs  where  L  =  20;  it  would  be 
identical  if  the  short  tank  length  were  L  =  12. 

In  conducting  these  simulations  using  shorter  tank  lengths,  it  is  discovered  that  a 
problem  arises  as  the  crest  approaches  the  end  of  the  tank.  The  original  Lagrangian  points 
closest  to  the  wall  tend  to  be  convected  away  from  the  wall  so  that  resolution  of  the  free 
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surface  near  the  wall  and  in  the  lower  part  of  the  wave  curl  becomes  poor.  In  order  to 
improve  the  resolution  in  this  regime,  a  routine  is  developed  to  add  points  in  the  wave  curl 
region  based  on  a  quadratic  interpolation  scheme  similar  to  that  used  in  the  original 
regridding  routine  for  the  free  surface  points.  This  technique  is  found  to  introduce  new 
difficulties  with  excessive  upward  motion  of  the  wall  contact  point  as  adjacent  points  are 
added.  Another  problem  arises  because  the  increased  density  of  Lagrangian  points  upon 
point  addition  in  some  cases  allows  excessive  "clustering"  of  points  in  the  curl  as  the 
simulation  proceeds  since  the  points  still  tend  to  be  convected  away  from  the  wall  and  into 
into  the  center  of  the  curl.  The  approaches  used  to  eliminate  or  alleviate  these  difficulties 
are  also  addressed. 

5.1  Results:  L=I  1.667 

5.1.1  Summary  of  Attempted  Runs  and  Failure  Types 

The  plots  of  free  surface  profiles  for  the  runs  referred  to  below  are  located  in 
Appendix  B. 

Run  B1  is  made  using  regridding  alone  every  15  time  steps.  Like  run  A1  in  Chapter 
4,  the  simulation  breaks  down  at  the  first  peak  when  the  overturning  process  is  not 
adequately  resolved.  The  time  of  failure  is  later  at  t  =  47.68  (2736)  around  x  =  1 1.0. 

Run  B2  is  made  using  smoothing  alone  every  15  time  steps  as  a  comparison  with  run 
Bl.  This  run  fails  in  a  manner  similar  to  that  in  run  A2  of  Chapter  4  when  the  first  peak 
begins  to  overturn.  A  thin  cusp  is  formed  which  eventually  self- intersects  before  re-entry, 
causing  failure  at  t  =  43.97  (5769)  around  x  -  9.3.  This  weak  breaking  event  at  the  first 
peak,  as  in  Chapter  4,  necessitates  some  means  of  artificial  suppression  in  order  to  allow  the 
second  peak  to  develop  to  the  point  of  breaking  near  the  end  of  the  tank. 

Run  B3  is  a  restart  of  run  B  l  at  t  =  42.00  (1800)  with  smoothing  every  5  time  steps. 
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It  is  analagous  to  run  A2  of  Chapter  4.  Failure  occurs  at  t  -  43.50  (3865)  at  the  first  peak 
when  the  overturning  crest  re-enters  the  free  surface  at  about  x  =  9.0.  Qualitatively,  runs 
B2  and  B3  are  similar  in  the  formation  of  a  thin  cu^p.  The  overturning  in  run  B3  is  delayed 
due  to  the  initial  phase  of  regridding. 

Run  B4  starts  the  simulation  from  t  =  0,  but  with  no  smoothing  or  regridding  used. 
The  RK4  scheme  followed  by  the  ABM  scheme  is  repeated  every  15  time  steps  (M  =  15) 
when  the  time  step  control  routine  is  called  to  compute  a  new  A /.  This  ran  is  of  particular 
interest,  as  its  objective  is  to  confirm  the  development  of  a  "sawtooth"  type  of  instability  in 
the  absence  of  any  smoothing  or  regridding.  This  type  of  instability  does  indeed  occur  in 
the  vicinity  of  the  wavemaker  and  causes  failure  at  t  =  31.7  (1476).  The  profiles  show  the 
development  of  the  instability,  which  is  first  detectable  at  the  scale  of  the  plots  at  about  t  = 
27.0(1000). 

Run  B5  is  a  restart  of  run  B1  at  t  -  42.00  (1800)  without  using  regridding  or 
smoothing.  As  in  ran  B4  the  time  integration  returns  to  the  RK4  scheme  every  15  time 
steps  when  the  time  step  size  is  recomputed.  Although  no  sawtooth  instability  is  observed 
in  this  run,  failure  occurs  at  the  first  peak  at  t  -  *3.88  (4639)  ac  the  crest  overturns  and  the 
cusp  self- intersects  before  touching  the  free  surface.  This  run  restarts  at  the  same  point  as 
run  B3,  and  a  comparison  shows  that  the  cusp  is  much  sharper  and  less  physical  in 
appearance  in  the  case  without  smoothing  applied. 

Run  B6  restarts  run  B1  at  r  =  42.00  (1800)  and  proceeds  Using  regridding  plus 
smoothing  every  5  time  steps  in  an  attempt  to  suppress  the  overturning  tendency  of  the  first 
peak  observed  in  runs  B2  through  B5  above.  It  is  analogous  to  run  A3  in  Chapter  4,  and  the 
result  is  the  same  in  that  overturning  of  both  the  first  and  second  peaks  is  suppressed.  This 
run  does  eventually  break  down  at  about  r  =  52.01  (3380)  as  the  second  peak  "sloshes"  up 
the  right  wall. 
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Runs  B6A  through  B6F  are  restarts  of  run  B6  at  the  following  respective  times  with 
smoothing  alone  every  5  steps:  43.84  (2000),  45.03  (2130),  45.99  (2215),  46.94  (2280)  and 
49.01  (2420).  As  in  runs  A3A  -  A3F  and  A4A  -  A4D  in  Chapter  4,  the  objective  in  these 
runs  is  to  capture  the  overturning  of  the  second  peak  after  having  suppressed  the 
overturning  of  the  Fust  peak  through  the  use  of  combined  regridding  and  smoothing  every  5 
steps. 

Run  B6A  fails  at  t  =  44.76  (2840)  when  the  overturning  second  peak  self-intersects 
before  touching  the  free  surface.  It  is  important  to  note  that  the  overturning  occurs  around  x 
=  9.8,  well  before  the  end  of  the  tank.  It  appears  as  though  there  is  insufficient  resolution  in 
the  curl  of  the  wave,  which  may  cause  errors  to  grow  in  this  region  of  high  velocities  and 
velocity  gradients  and  lead  to  the  observed  failure  at  the  top  of  the  curl. 

Run  B6B  fails  at  t  =  49.02  (2922)  when  the  overturning  peak  self-intersects  before 
touching  the  free  surface.  In  this  case  also  it  appears  that  there  is  insufficient  density  of 
Lagrangian  points  along  the  crest  tip  and  curl  to  accurately  capture  the  overturning.  The 
overturning  occurs  later  than  run  B6A  and  closer  to  the  right  wall,  around  x  =  10.4. 

Run  B6C  fails  at  t  =  5 1 .06  (3740)  when  the  overturning  crest  again  self-intersects  in  a 
manner  similar  to  runs  B6A  and  B6B.  The  overturning  occurs  even  closer  to  the  right  wall, 
around  x  =  11.63. 

Run  B6D  is  identical  to  B6C  except  that  the  minimum  allowable  time  step  size  is  set 

\ 

to  0.00001  instead  of  0.0001 ,  which  was  used  in  all  of  the  above  cases.  The  purpose  of  this 
is  to  see  if,  in  the  overturning  phase,  reducing  the  restriction  on  time  step  size  will  affect  the 
results.  Although  more  time  steps  are  required  because  of  reduced  Ar,  failure  still  occurs  in 
essentially  the  same  fashion  at  t  =  5 1 .06  (3766).  No  profiles  are  plotted  for  this  case. 

Run  B6E  continues  until  the  overturning  crest  is  almost  at  the  point  of  re-entry  when 
failure  occurs.  One  difference  between  this  run  and  the  previous  runs  is  that  the  density  of 
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Lagrangian  points  at  the  crest  tip  and  in  the  curl  is  somewhat  higher,  and  could  explain  the 
fact  that  the  run  does  not  fail  as  soon;  however,  the  curl  does  distort  in  a  way  similar  to  runs 
6A  -  6D  just  prior  to  failure.  The  crest  tip  is  just  above  the  free  surface  at  approximately  t  = 
51.034  (3792)  and  x  =  11.63.  The  panel  adjacent  to  the  wall  is  0.7  times  the  size  of  the 
original  panel  size  at  t  -  0  by  this  point.  The  panels  decrease  in  size  with  distance  from  the 
right  wall,  and  reach  a  minimum  size  at  the  crest  tip  of  0.08  times  the  original  panel  size. 
Failure  occurs  at  t-  51.035  (3802). 

Run  B6F  fails  at  t  =  50.96  (3662)  as  the  overturning  crest  self-intersects  around  x  = 
11.56.  Compared  to  run  B6E,  the  resolution  of  the  curl  region  is  worse  by  roughly  one 
panel,  and  distortion  at  the  top  of  the  curl  similar  to  runs  B6A  -  B6D  continues  until  failure 
occurs. 

Run  B7  duplicates  run  B6  except  that  smoothing  is  applied  just  before  regridding 
rather  than  just  after.  The  switch  has  no  appreciable  effect  on  the  results,  and  profile  plots 
are  not  given. 

5.1.2  Comparison  of  Successful  Run  with  Experiment 

From  the  above  runs,  B6E  is  considered  successful  since  the  simulation  fails  just 
before  re-entry  occurs,  even  though  resolution  of  the  wave  curl  region  is  questionable.  As 
anticipated,  overturning  does  occur  earlier  than  in  the  actual  experiment,  and  the  crest  tip 
nears  re-entry  before  the  crest  reaches  the  right  wall.  At  this  point  it  is  of  interest  to 
compare  free  surface  elevations  with  the  profiles  given  for  the  actual  experiment  prior  to 
impact  by  Chan  [2]  in  Figure  6  for  the  wall  position  (tank  length)  of  L  =  11.667.  At  t 
=51.49,  the  maximum  crest  height  is  about  0.18  and  the  run-up  at  the  right  wall  is  about 
0.13  just  as  overturning  begins  At  t  =  51.53,  just  before  impact  occurs,  the  maximum  crest 
height  is  about  0.21  and  the  run-up  is  about  0.16.  In  the  simulation,  just  as  overturning 
begins,  the  maximum  crest  height  is  about  0.17  and  the  wall  run-up  is  about  0.10.  The 
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maximum  crest  height  at  the  end  of  the  non  is  about  0.17  and  the  run-up  at  the  right  wall  is 
about  0.11.  The  crest  height  and  run-up  for  the  simulation  are  thus  comparable  to 
experiment  as  overturning  starts,  but  because  the  simulated  wave  never  reaches  the  wall  the 
rapid  increase  in  crest  height  and  elevation  at  the  wall  observed  in  the  experiment  just  prior 
to  impact  does  not  occur. 


In  the  simulation,  the  height  of  the  curl  region  as  the  wave  overturns  is  on  the  order 
of  0.01  -  0.02.  Assuming  that  a  simulation  can  eventually  be  achieved  where  overturning 
has  started  but  re-entry  has  not  occurred  as  the  crest  reaches  the  wall  boundary,  the  height 
of  the  space  "trapped"  between  the  impacting  crest  and  free  surface  at  the  wall  at  the  bottom 
of  the  curl  should  be  cr>  the  same  order.  Using  the  same  profiles  from  the  actual  experiment 
an  estimate  of  the  trapped  air  pocket  height  at  impact  is  about  0.02  -  .04  for  L  =  1 1.667. 


In  order  to  achieve  the  desired  simulated  result,  an  obvious  choice  is  to  try  decreasing 
the  tank  length  until  the  overturning  wave  reaches  the  wall  before  re-entry.  This  approach 
is  described  in  the  next  section. 


5.2  Results:  Variation  of  Tank  Lengths 

In  the  following  runs,  the  same  sequence  of  regridding  and  smoothing  is  used  as  for 
run  B6E.  The  initial  run  is  made  using  regridding  alone  every  15  steps  until  Trs  =  42.0, 
followed  by  regridding  plus  smoothing  every  5  steps  until  Ts  =  47.0,  after  which  smoothing 
alone  is  used  every  5  steps  until  the  end  of  the  simulation.  In  all  runs  the  minimum 
allowable  time  step  size  is  set  to  0  0001.  Plots  of  the  final  free  surface  profiles  for  these 
runs  are  located  in  Appendix  C 


t 
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5.2.1  Summary  of  Attempted  Runs  and  Failure  Types 


The  following  table  summarizes  the  runs  for  varying  tank  length. 


TABLE  5.1 

RUN 

LENGTH 

Nrs 

Ns 

Nf 

Tf 

Cl 

usm 

1815 

2310 

3282 

ESS 

C2 

1815 

2305 

3350 

EEJ 

C3 

11.59 

1810 

2300 

3525 

50.95 

C4 

11.60 

1800 

2295 

3590 

50.97 

C5 

11.63 

1805 

2290 

3652 

51.00 

C6 

11.65 

1800 

2285 

3740 

50.02 

Nrs  and  Ns  axe  the  time  steps  corresponding  to  Ts  and  Trs;  Tf  and  Nf  are  the  time  and  time 
step  at  which  the  run  fails. 

Runs  Cl  and  C2  fail  in  a  similar  fashion.  In  these  cases,  the  tank  length  is  so  short 
that  the  free  surface  at  the  wall  rises  rapidly  relative  to  the  forward  motion  of  the  wave 
crest;  as  a  result,  overturning  never  occurs.  The  runs  fail  when  the  point  adjacent  to  the 
wall  contact  point  is  convected  to  the  left  and  eventually  overtakes  the  points  to  its  left. 


In  runs  C3  and  C4,  overturning  does  begin  due  to  longer  tank  lengths.  The  profiles 
for  run  C3  are  very  similar  to  run  C4  at  the  start  of  overturning,  only  earlier  in  time  by 
about  0.016.  At  comparable  points  of  t  =  50.94  (3425)  for  run  C3  and  t  =  50.96  (3500)  for 
run  C4,  just  after  overturning  has  begun,  the  panel  next  to  the  right  wall  is  about  0.4  times 
the  original  panel  length.  Both  simulations  eventually  fail  as  the  point  adjacent  to  the  wall 
point  moves  so  far  to  the  left  that  there  is  no  resolution  at  all  of  the  the  lower  part  of  the 
curl.  Accurate  simulation  of  overturning  is  no  longer  possible  and  the  curl  is  distorted  in 
such  a  way  that  the  crest  "folds"  down  on  itself  instead  of  plunging. 


In  runs  C5  and  C6  overturning  also  begins,  but  in  these  cases  the  tank  is  too  long 
because  the  crest  is  too  far  from  the  wall.  They  fail  in  a  way  similar  to  runs  B6A  -  B6D  and 
run  B6F  in  section  5.1.1.  The  overturning  in  these  cases  is  later  than  in  runs  C3  and  C4, 
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and  it  is  interesting  to  note  that  in  these  four  cases  the  overturning  time  is  delayed  as  tank 
length  increases. 

Despite  the  failure  of  runs  C3  and  C4,  they  are  still  potentially  useful.  Before  failure, 
the  trajectory  of  the  overturning  crest  tip  is  roughly  horizontal  and  close  enough  to  the  wall 
to  suggest  that  solving  the  problem  with  resolution  near  the  wall  in  these  cases  may  allow 
the  simulation  to  continue  to  impact.  Run  C4  with  L  =  11.60  is  chosen  to  continue  efforts 
to  achieve  an  impact  case. 

5.2.2  Comparison  of  Run  L=1 1.60  with  Experiment 

At  r  =  50.958  (3500),  the  profile  of  run  C4  shows  a  maximum  crest  height  of  0.18 
and  wall  run-up  of  0. 15.  This  compares  reasonably  well  to  the  profile  of  Chan  [2]  for  L  = 
1 1.667  and  t  =  51.49  which  shows  a  maximum  crest  height  of  about  0.18  and  wall  run-up  of 
about  0.13  at  roughly  the  same  stage  of  overturning.  Before  failure  occurs  it  appears  that 

i 

the  simulation  might  continue  to  form  a  "pocket"  between  the  overturning  crest  and  wall  at 
impact  with  an  initial  height  of  about  0.01. 

1  5.3  Results:  Effort  to  Increase  Point  Resolution 

5.3.1  Techniques  Investigated 

i  The  following  sections  describe  the  approach  that  is  taken  to  address  the  problem  of 

insufficient  resolution  near  the  wall  and  in  the  curl.  Run  C4  is  used  as  the  baseline  run 
from  which  changes  are  made. 

1  5.3. 1.1  Point  addition 

Addition  of  Lagrangian  points  is  the  first  technique  used  to  increase  point  resolution 
near  the  wall  and  in  the  wave  curl.  A  subroutine  is  developed  which  first  finds  the  smallest 
'  panel  on  the  free  surface  and  then  adds  points  at  the  midpoints  of  subsequent  panels  when 
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their  length  exceeds  a  specified  multiple  of  the  minimum  panel  length;  in  these  runs  two  is 
the  multiple.  The  points  are  added  using  quadratic  interpolation  for  position  and  velocity 
potential  with  the  linear  distance  between  panels  as  the  parameter.  This  subroutine  is  called 
immediately  before  the  smoothing  subroutine  is  called,  in  these  cases  every  five  time  steps. 

In  the  first  series  of  runs,  run  C4  is  restarted  at  t  =  50.958  (3500)  just  after 
overturning  has  begun.  The  free  surface  profiles  for  these  runs  are  located  in  Appendix  D. 

In  run  Dl,  smoothing  is  used  every  5  steps  and  points  addition  is  applied  once, 
causing  the  addition  of  3  points  at  step  3505  at  the  midpoints  of  the  three  panels  which  are 
nearest  the  right  wall.  The  miminum  panel  size  used  for  comparison  in  the  point  addition 
routine  occurs  at  the  tip  of  the  overturning  crest  where  the  Lagrangian  points  have 
concentrated,  and  is  0.05  times  the  original  panel  size. 

Examining  the  sequence  of  profiles  for  run  Dl  from  steps  3510  to  3570,  and 
comparing  them  with  the  same  sequence  of  profiles  for  run  C4  without  point  addition,  it  is 
apparent  that  point  addition  has  improved  resolution  in  part  of  the  curl  farthest  from  the 
wall  but  at  the  same  time  caused  the  contact  point  at  the  right  wall  to  move  rapidly  upward 
in  a  non-physical  way  with  an  average  velocity  of  about  4.0  until  the  overturning  crest  tip 
contacts  this  upward-moving  surface  and  the  run  fails. 

Run  D2  is  identical  to  Dl  except  that  the  minimum  allowable  time  step  size  is  set  to 
0.00001  instead  of  0.0001,  because  in  run  Dl  the  miminum  time  step  size  was  reached. 
The  new  minimum  time  step  is  never  reached  but  the  results  of  the  run  are  indistiguishable 
from  run  Dl  at  a  given  time  value,  eliminating  time  step  size  restriction  as  a  possible  cause 
of  the  problem. 

Run  D3  is  an  extreme  case  in  that  point  addition  is  applied  10  times  starting  at  step 
3505  and  ending  at  3600,  resulting  in  a  total  addition  of  15  points  in  the  curl-wall  region. 
Examining  the  sequence  of  plots  from  step  3515  to  3575  shows  that  the  rapid  upward 


-45- 


movement  at  the  wall  is  more  pronounced  and  more  localized  near  the  wall  with  the 
addition  of  points,  and  the  run  fails  sooner  than  runs  Dl  and  D2  due  to  the  anomalous 
behavior  of  the  points  nearest  the  wall  in  the  "jet". 

Run  D4  is  identical  to  runs  Dl  and  D2  except  that  no  smoothing  is  applied  after 
restart.  It  is  important  to  note  here  that  the  smoothing  routine  used  thus  far  smooths  not 
only  interior  points  on  the  free  surface  but  the  contact  points  as  well  using  a  one-sided 
formula  instead  of  a  centered  formula.  The  values  smoothed  are  the  position  and  velocity 
potential.  Comparing  the  profile  of  run  Dl  at  step  3570  with  run  D4  at  3565,  the  position 
of  the  contact  point  is  almost  identical.  The  crest  tip  in  run  D4  is  thinner  and  turned  more 
downward  than  in  Dl,  and  there  is  evidence  of  instability  in  run  D4  due  to  lack  of 
smoothing  in  the  slightly  jagged  appearance  of  points  at  the  crest  tip  and  in  the  curl.  For 
the  above  cases,  the  conclusion  is  that  a  jet  at  the  wall  is  formed  regardless  of  whether 
smoothing  is  applied. 

The  following  five  runs  restart  run  C4  earlier,  at  t  -  50.846  (3100)  before  overturning 
occurs  to  see  if  there  is  any  dependence  on  the  time  of  point  addition. 

Run  D5  uses  smoothing  every  5  steps  after  restart  and  applies  point  addition  once, 
resulting  in  the  addition  of  seven  points  at  step  3 105.  Examining  the  profiles  for  steps  3130 
to  3530  reveals  that  although  there  is  still  a  local  upward  deformation  of  the  free  surface  at 
the  wall  after  point  addition,  it  is  less  pronounced  with  an  average  velocity  of  about  1.0 
over  this  range.  As  the  wave  overturns  and  approaches  the  wall,  the  disparity  in  elevation 
of  the  contact  point  and  adjacent  points  becomes  less  marked.  From  steps  3630  to  3930,  (r 
=  50.960  to  t-  50.973)  the  overturning  process  appears  much  more  physical  than  in  run  C4 
over  a  similar  range  from  steps  3500  to  3580  (r  =  50.958  to  t  =  50.966)  with  no  points 
added.  However,  failure  of  run  D5  eventually  does  occur  in  a  way  similar  to  run  C4  when 
the  point  next  to  the  wall  contact  point  begins  to  move  to  the  left  and  eventually  leaves  no 
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resolution  of  the  region  near  the  wall;  the  crest  then  "folds"  down  on  itself  as  the  crest  tip 
reaches  the  wall.  The  time  of  failure  is  later  than  run  C4  by  about  0.006. 

Run  D6  is  identical  to  run  D5  except  point  addition  is  applied  twice,  resulting  in 
addition  of  12  points,  seven  at  step  3105  and  five  at  3110.  The  density  of  points  near  the 
wall  is  doubled  relative  to  run  D5  and  quadrupled  relative  to  run  C4.  Examining  the 
profiles  from  steps  3145  to  3445  reveals  that  the  upward  deformation  at  the  wall  has  an 
initial  average  velocity  of  about  l  .0,  as  in  run  D5.  This  run  also  fails  because  of  eventual 
interference  of  the  point  next  to  the  wall  with  other  points  as  it  moves  to  the  left,  causing 
loss  of  resolution  near  the  wall  and  deformation  at  the  bottom  of  the  curl.  The  time  of 
failure  is  about  0.004  later  than  run  C4. 

Run  D7  is  identical  to  run  D6  except  point  addition  is  applied  three  times,  resulting  in 
the  addition  one  more  point  at  step  3255  for  a  total  of  13  points  added  instead  of  12.  The 
results  are  very  similar  to  run  D6. 

The  next  two  runs  are  made  to  see  what  affect  no  smoothing  has  on  the  results. 

Run  D8  is  identical  to  run  D5  except  that  smoothing  is  not  used;  point  addition  is 
applied  once,  resulting  in  the  addition  of  seven  points.  Examining  the  profiles  for  steps 
3101  to  3501,  there  is  still  an  upward  deformation  of  the  free  surface  at  the  wall  which 
again  has  an  average  velocity  of  about  1.0  over  this  range.  The  striking  aspect  of  this  ran, 
however,  is  the  failure  due  to  the  sawtooth  instability  which  develops  in  the  curl  region  and 

i 

near  the  wall  as  the  wave  overturns.  In  this  case  failure  occurs  just  below  the  crest  tip  at 
one  of  the  original  points  having  an  added  point  on  either  side.  It  is  not  possible  from  this 
ran  alone  to  say  what  affect  point  addition  has  on  the  development  of  the  sawtooth 
instability.  The  time  of  failure  is  about  0.018  sooner  than  run  C4. 

Run  D9  is  identical  to  run  D7  except  that  smoothing  is  not  used;  point  addition  is 
applied  three  times,  resulting  in  the  addition  of  13  points.  Sawtooth  instability  is  also 
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evident  in  the  wave  curl,  but  this  run  fails  about  0.01  earlier  than  run  D4  due  to  an  anomaly 
at  the  wall  where  the  the  characteristic  jet  has  developed.  The  time  of  failure  is  about  0.026 
sooner  than  run  C4. 

In  summary,  runs  D1  through  D9  reveal  the  following: 

1 .  Point  addition  causes  the  velocity  of  the  wall  contact  point  to  increase  in  a 
non-physical  way,  causing  upward  deformation  of  the  free  surface  at  the  wall. 

2.  The  formation  of  this  "jet"  is  dependent  on  when  point  addition  occurs.  In 
runs  D1  -  D4  the  distortion  of  the  jet  is  most  severe  and  clearly  non-  physical; 
in  runs  D5  -  D9  the  distortion  is  less  pronounced  when  points  are  added  near 
the  wall  point  earlier  in  the  overturning  process. 

3.  This  jet  is  formed  whether  smoothing  is  applied  or  not. 

4.  In  cases  D4  -  D7,  where  the  jet  at  the  wall  is  less  dominant  an  effect,  the 
simulation  is  improved  relative  to  run  C4  due  the  better  resolution  in  the  curl; 
however,  as  in  run  C4,  leftward  movement  of  the  point  next  to  the  wall  point 
eventually  causes  loss  of  resolution  near  the  wall  and  distortion  of  the  curl, 
leading  to  failure. 

5.  Because  of  relatively  early  failure  of  runs  D8  and  D9,  it  is  not  yet  clear  if  the 
motion  of  this  point  to  the  left  is  affected  by  the  smoothing  of  points  near  the 
wall. 

6.  Sawtooth  instabilities  do  occur  if  no  smoothing  is  applied  upon  restart. 

The  difficulty  with  the  leftward  movement  of  the  point  closest  to  the  wall  is  apparent 
both  before  and  after  point  addition  is  attempted,  whereas  the  formation  of  the  jet  at  the 
wall  occurs  only  after  point  addition. 
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5.3.1.2  Modified  smoothing 

Because  both  problems  involve  the  motion  of  the  wall  contact  point  and  immediately 
adjacent  points,  it  is  conceivable  that  the  application  of  smoothing  may  need  to  be  modified 
in  some  way.  As  demonstrated  in  runs  D8  and  D9,  some  amount  of  smoothing  is  clearly 
necessary  in  order  to  prevent  sawtooth  instability.  The  original  smoothing  routine  smooths 
the  position  and  velocity  potential  of  all  points  on  the  free  surface,  including  the  wall 
contact  points.  For  the  interior  points,  a  centered  5-point  formula  is  used;  for  the  point 
adjacent  to  the  wall  point,  a  skewed  5-point  formula  is  used;  and  for  the  wall  point  a  one¬ 
sided  5-point  formula  is  used.  The  smoothed  values  at  each  point  are  computed  using  the 
original  (unsmoothed)  values  at  adjacent  points,  and  subsitution  of  the  new  values  for 
original  values  is  done  as  the  final  step. 

The  runs  in  this  section  investigate  the  effect  of  variations  in  use  of  the  original 
smooting  routine  and  changes  to  the  original  routine  with  regard  to  how  the  wall  points  and 
adjacent  points  are  handled.  The  free  surface  profiles  for  these  runs  are  located  in 
Appendix  E. 

All  these  runs  have  identical  input  to  run  D6:  they  restart  ran  C4  at  step  3500  and 
apply  point  addition  twice,  adding  seven  points  at  step  3505  and  five  points  at  3510. 

In  run  El,  the  order  of  subroutines  is  switched  so  that  smoothing  of  free  surface 
points  occurs  first  followed  by  the  addition  of  any  points.  The  affect  upon  the  free  surface 
profiles  is  negligible  when  compared  to  run  D6.  Plots  of  profiles  are  not  given  for  this  run. 

In  run  E2,  the  original  smoothing  routine  is  used  but  the  arguments  are  changed  so 
that  the  values  of  velocity  potential  and  position  at  the  wall  point  are  neither  smoothed 
themselves  nor  used  in  the  formula  for  smoothing  the  adjacent  points;  in  other  words,  the 
point  adjacent  to  the  wall  is  taken  as  the  end  point  and  the  wall  point  is  ignored.  The  results 
of  this  run  are  significantly  different  from  run  D6.  In  examining  the  profiles  from  step 
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3161  to  3861,  note  that  the  wall  point  still  moves  rapidly  upward  after  point  addition,  but 
that  the  adjacent  points  do  not.  By  ignoring  the  wall  point  in  using  the  smoothing  routine, 
the  motion  of  the  adjacent  points  is  no  longer  coupled  "artificially"  to  the  wall  point  motion 
through  smoothing,  though  there  is  still  a  coupling  to  the  wall  point  motion  in  the  solution 
of  the  boundary  value  problem  and  subsequent  computation  of  velocities.  The  adjacent 
points  do  eventually  "catch  up"  to  the  wall  point,  as  seen  at  step  3761  when  the  crest  tip  is 
about  0.01  from  the  wall,  then  the  wall  point  again  moves  upward  relative  to  the  adjacent 
points. 

Another  significant  difference  from  run  D6  is  seen  in  the  behavior  of  the  point 
adjacent  to  the  wall.  In  run  D6  the  movement  of  this  point  to  the  left  becomes  noticable 
around  step  3445  (r  =  50.946),  and  it  eventually  interferes  with  the  points  to  its  left  to  cause 
deformation  of  the  bottom  of  the  curl  and  failure.  By  step  3645  (r  =  50.960),  for  instance, 
the  adjacent  point  to  the  wall  is  extremely  close  to  the  second  point  from  the  wall, 
producing  a  panel  of  length  0.001,  about  37  times  smaller  than  the  original  panel  size,  while 
the  panel  between  the  wall  and  the  adjacent  point  has  a  length  of  0.128,  about  0.3  times  the 
size  of  the  original  panel  (0.039).  In  contrast,  tun  D6  at  step  3661  (r  =  50.961)  shows  no 
large  deviation  in  position  of  the  point  next  to  the  wall.  The  panel  next  to  the  wall  has 
length  0.018  and  the  second  panel  from  the  wall  has  length  0.015.  Although  the  point 
closest  to  the  wall  is  moving  to  the  left,  the  motion  is  not  nearly  as  apparent. 

The  failure  of  run  E2  is  characterized  by  the  clustering  together  of  points  into  the 
center  of  the  curl  until  the  spacing  becomes  extremely  small  and  computation  breaks  down. 
Comparison  of  the  profiles  at  steps  3861  and  3961  shows  the  motion  of  points  at  the  bottom 
of  the  curl  into  the  center.  Note  also  that  the  wall  point  is  again  significantly  higher  than  the 
adjacent  point  by  the  end  of  the  run,  as  it  was  at  the  beginning,  and  that  the  adjacent  point 
has  moved  noticeably  away  from  the  wall  but  not  so  far  as  to  cause  the  kind  of  failure  in 
run  D6.  Failure  occurs  at  t  =  50.975,  about  0.004  later  than  run  D6. 
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In  summary,  the  omission  of  the  wail  point  in  use  of  the  original  smoothing  routine 
leads  to  the  following  observations: 

1.  The  increase  in  wall  point  velocity  upon  point  addition  still  occurs,  but  unlike 
previous  runs  the  adjacent  points  do  not  move  upward  with  it.  The  wall  point 
motion  seems  effectively  independent  of  the  adjacent  point  motion  when  the 
wall  point  is  not  used  as  the  end  point  in  the  smoothing  computation. 

2.  The  simulation  proceeds  considerably  farther  than  ran  D6  due  to  the  fact  ihat 
the  point  adjacent  to  the  wall  does  not  move  rapidly  to  the  left.  The 
overturning  crest  tip  even  passes  the  end  of  the  tank  since  the  program  has  not 
yet  been  modified  to  define  impact  at  the  end  of  the  tank. 

3.  The  addition  of  points  to  increase  the  resolution  in  the  curl  and  near  the  wall 
eventually  causes  difficulty  as  the  simulation  proceeds  because  of  the 
observed  tendency  of  the  points  at  the  bottom  of  the  curl  and  near  the  wall  to 
move  toward  the  center  of  the  curl,  causing  failure  when  the  point  spacing 
becomes  extremely  small. 

The  fact  that  the  anomalous  upward  acceleration  of  the  wall  point  occurs  whether  or 
not  the  wall  point  is  used  in  smoothing,  along  with  the  fact  that  the  adjacent  point  moves 
rapidly  to  the  left  only  when  smoothing  does  include  the  wall  point,  suggest  that  non¬ 
physical  motion  of  the  wall  point  might  induce  error  in  the  adjacent  point  velocity  when 
coupled  through  the  smoothing  operation. 

In  run  E3,  the  original  smoothing  routine  is  again  used.  In  this  case,  the  position  of 
all  points  is  smoothed  (except  the  x-position  of  the  wall  point,  which  can  only  move 
vertically)  but  the  velocity  potential  is  not  smoothed.  For  this  run,  each  plot  of  the  free 
surface  profile  is  followed  by  a  corresponding  plot  of  velocity  potential  versus  x-position. 
Examination  of  plots  for  steps  3150  to  3230  show  that  a  jet  is  again  formed  by  the  wall 
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point  and  adjacent  points,  but  that  upward  deformation  is  greater  than  in  runs  where 
velocity  potential  end  position  were  smoothed.  At  step  3230  (/  =  50.899),  just  as  the  crest 
is  starting  to  overturn,  note  that  the  wall  point  elevation  is  comparable  to  the  crest  height 
and  that  the  tip  of  the  crest  is  not  "smooth"  but  a  sharp  comer.  There  is  a  corresponding 
distortion  in  the  velocity  potential  plot  at  the  crest  tip  and  wall.  The  free  surface  points  also 
have  a  slightly  jagged  appearance,  suggesting  that  instability  may  result  when  smoothing  of 
velocity  potential  does  not  occur,  even  if  smoothing  of  position  does  occur.  This  run 
eventually  fails  at  t  =  50.907,  sooner  than  run  D6  by  about  0.065. 

In  run  E4,  the  smoothing  routine  itself  is  modified  so  that  the  wall  point  position  and 
velocity  potential  are  not  smoothed  themselves  but  are  used  in  the  smoothing  formula  for 
the  adjacent  two  points  This  run  differs  from  run  E2,  where  the  wall  point  values  are  not 
used  at  all  in  smoothing.  Examining  profiles  for  steps  3183  to  3983,  the  initial  jet  at  the 
wall  again  develops,  but  another  problem  is  observed  to  occur  when  the  fourth,  fifth  and 
sixth  points  from  the  wall  become  very  closely  spaced  as  the  points  closest  to  the  wall  move 
to  the  left.  This  leads  to  local  deformation  of  the  free  surface  at  these  points  and  eventual 
failure  around  t  =  50.968.  The  degree  of  deformation  can  be  seen,  for  instance,  by 
comparing  the  profile  at  step  3983  (t  =  50.967)  with  run  E2  at  step  3761  (t  =  50.967). 

Although  the  jet  at  the  wall  still  occurs  in  run  E4,  it  is  significant  to  note  that  the 
point  closest  to  the  wall  does  not  move  as  rapidly  to  the  left  as  it  does  in  run  D6,  but  that  the 
motion  still  is  greater  than  in  run  E2  where  the  wall  point  is  not  used  at  all  in  smoothing. 
This  observation  supports  the  conjecture  that  the  anomalous  wall  point  behavior  may 
corrupt  the  computation  of  values  at  the  adjacent  points  since  the  smoothing  operation 
couples  the  values  of  these  points,  even  though  in  run  E4  the  wall  point  values  themselves 
are  not  smoothed.  Conversely,  the  adjacent  point  motion  could  significantly  affect  the  wall 
point  motion  through  the  smoothing  operation. 
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5.3. 1.3  Alternate  finite  difference  for  wall  point  velocity 

The  increase  in  wall  point  velocity  upon  point  addition,  evident  from  the  plots  of  free 
surface  profile,  is  confirmed  by  checking  actual  velocity  values.  For  example,  in  run  C4, 
which  does  not  add  points,  the  change  in  vertical  velocity  at  the  wall  point  from  step  3104 
to  3105  is  0.003,  and  from  step  3105  to  step  3106  is  0.002.  In  run  D6,  which  adds  seven 
points  at  step  3105  then  smooths  values,  the  respective  changes  in  wall  point  velocity  are 
0.064  and  0.003.  In  run  D9,  which  adds  seven  points  at  3105  but  does  not  smooth  any 
values,  the  respective  changes  in  wall  point  velocity  are  0.062  and  0.003.  In  run  E2,  which 
excludes  the  wall  point  from  smoothing,  seven  points  are  also  added  at  3105,  and  the 
respective  changes  in  wall  point  velocity  are  0.064  and  0.003.  In  all  cases  of  point  addition 
the  jump  in  velocity  from  step  3104  to  3105  is  about  eight  per  cent,  compared  to  a  jump  of 
0.4  per  cent  in  run  C4  without  point  addition.  From  step  3105  to  3106  the  change  in 
velocity  returns  to  a  value  comparable  to  run  C4.  When  points  are  added  again  at  step 
3110,  there  is  again  a  jump  in  velocity  from  step  3109  to  3110  of  about  sixteen  percent 
based  on  the  average  velocity.  For  the  range  of  time  steps  under  consideration  the  time  step 
size  is  about  0.0007. 

Note  that  in  the  original  and  modified  programs  used  thus  far,  the  wall  contact  point 
velocity  is  computed  by  a  one-sided  three-point  difference  formula  applied  to  the  complex 
potential  along  the  free  surface.  Therefore,  when  a  point  is  added  between  the  wall  point 
and  the  original  point  adjacent  to  the  wall,  the  subsequent  computation  of  contact  point 
velocity  is  now  based  in  part  on  the  interpolated  position  and  velocity  potential  at  this  new 
point.  The  magnitude  of  the  difference  Az  used  in  the  difference  formula  is  also  cut  in  half 
upon  addition  of  the  point.  If  a  second  point  is  later  added  between  the  wall  point  and  the 
previous  new  point,  as  it  is  in  the  cases  above,  computation  of  the  wall  point  velocity  no 
longer  uses  any  original  adjacent  points,  and  the  magnitude  of  Az  is  again  cut  in  half.  It 
should  be  pointed  out  that  the  identities  of  the  points  used  in  the  smoothing  routines 
discussed  so  far  also  change  when  points  are  added. 
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It  is  possible  that  point  addition  could  cause  the  observed  jumps  or  errors  in  the  wall 
point  velocity  because  of  the  introduction  of  quadratically  interpolated  values  and  smaller 
point  spacing  into  the  present  formula  for  computing  wall  point  velocity.  The  computation 
of  velocity  at  interior  points  might  also  be  significantly  affected  when  new  adjacent  points 
are  added. 

It  is  interesting  to  recall  that  when  points  are  added  later  in  the  overturning  process 
the  deformation  at  the  wall  is  much  more  pronounced  that  when  points  are  added  earlier 
53.1.1.  This  suggests  that  any  error  introduced  in  the  wall  point  velocity  upon  point 
addition  is  larger  when  points  are  added  at  a  later  time,  corresponding  to  greater  velocities 
of  the  free  surface  near  the  wall. 

An  alternate  method  of  computing  the  vertical  velocity  of  the  wall  points  after  each 
solution  step  is  to  use  a  one-sided  finite  difference  of  the  velocity  potential  only  along  the 
right  wall,  since  d<| )/dy  =v.  This  approach  differs  significantly  in  that  the  points  used  are 
always  evenly  spaced,  points  are  never  added  to  this  region,  and  the  values  of  velocity 
potential  and  position  at  the  points  below  the  wall  point  are  not  directly  affected  by  the 
smoothing  of  the  free  surface  points. 

Based  on  this  approach,  3 -point,  4-point  and  5-point  finite  difference  formulas  are 
subsituted  for  the  original  scheme  when  computing  the  wall  point  vertical  velocity  in  the 
following  runs.  The  original  smoothing  routine  is  still  used  in  which  the  wall  point  values 
are  smoothed  as  well  as  the  interior  free  surface  points.  These  runs  restart  run  C4  at  t  = 
50.846  (3100)  and  apply  smoothing  every  five  steps.  All  the  runs  apply  point  addition 
twice,  adding  seven  points  at  step  3105  and  five  points  at  3110.  The  free  surface  profiles 
for  these  runs  are  found  in  Appendix  F. 

In  run  FI,  the  3-point  formula  is  used.  The  free  surface  profiles  show  no  significant 
change  in  the  vertical  motion  of  the  wall  point  compared  to  run  D6  which  uses  the  original 
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velocity  computation  at  the  wall  point.  This  run  tails  at  t  -  50.973  (3985)  in  a  way  very 
similar  to  run  D6  as  the  point  adjacent  to  the  wall  moves  rapidly  to  the  left  as  overturning 
proceeds,  eventually  interfering  with  adjacent  points  and  causing  deformation  at  the  bottom 
of  the  curl. 

Runs  F2  and  F3  use  the  4-point  and  5-point  formulas,  respectively.  Again,  there  is 
little  difference  in  the  outcome  compared  to  mn  D6.  The  failure  times  for  run  F2  and  F3 
are,  respectively,  t  -  50.973  (4046)  and  t  =  50.972  (4065).  No  profiles  are  plotted  for  these 
runs. 

5.3.1.4  Alternate  Finite  difference  plus  modified  smoothing 

The  new  method  of  computing  wall  point  velocity  appears  to  have  little  effect  on  the 
results  using  the  original  smoothing  technique.  The  following  runs  use  the  new  method  in 
combination  with  variations  of  smoothing  already  used  in  runs  E2  and  E4  with  the  original 
method  of  velocity  computation.  The  free  surface  profiles  for  all  runs  in  this  section  are 
found  in  Appendix  G. 

The  next  three  runs  change  the  smoothing  routine  to  prevent  smoothing  of  the  wall 
point  vertical  position  and  velocity  potential  but  do  use  these  values  in  smoothing  the 
adjacent  two  points.  This  smoothing  technique  is  identical  to  that  used  in  mn  E4,  in  which 
the  point  adjacent  to  the  wall  is  seen  to  move  left  much  less  than  in  mn  D6,  although  failure 
still  occurs  due  to  problems  with  the  clustering  of  points  near  the  wall.  These  mns  again 
restart  mn  C4  at  step  3100  and  apply  smoothing  every  five  steps.  Point  addition  is  applied 
twice,  adding  seven  points  at  step  3105  and  five  points  at  step  3110. 

In  mn  Gl,  which  uses  the  3-point  velocity  formula  for  the  wall  point,  there  is  a 
significant  difference  in  the  behavior  of  the  wall  point  and  adjacent  points.  Examining  the 
profiles  for  steps  31 10  to  4210  shows  that  initially  (i.e.,  through  step  3210)  the  wall  point 
and  adjacent  points  remain  essentially  even  in  elevation,  with  no  jet  formed.  From  steps 


3310  through  3810  the  wall  point  is  seen  to  lag  behind  the  adjacent  points  slightly  in 
upward  velocity  and  position.  Note  that  at  step  3810  the  crest  tip  has  nearly  reached  the 
boundary  of  the  wall,  and  the  simulation  looks  reasonable  except  for  the  low  postion  of  the 
wall  point  compared  to  the  adjacent  points.  It  is  also  significant  that  by  this  step  there  has 
been  no  marked  change  in  spacing  of  the  points  near  the  wall  which  characterize  the  later 
stages  of  previous  runs.  As  the  crest  tip  passes  the  wall  boundary,  an  anomaly  occurs  at  the 
wall  in  which  the  points  closest  to  the  wall  point  pass  over  the  wall  point,  creating  a 
deformation  which  eventually  leads  to  failure.  The  time  of  failure  is  t  =  50.995  (4410); 
however,  this  relatively  late  time  is  due  to  the  abnormal  evolution  of  the  free  surface  and  is 
not  a  good  basis  of  comparison  with  previous  runs. 

In  run  G2,  which  uses  the  4-point  formula,  the  profiles  for  steps  3152  through  4052 
show  that  the  difficulty  in  wall  point  motion  observed  in  run  G1  is  not  present.  The  wall 
point  never  exhibits  a  marked  difference  in  elevation  relative  to  the  adjacent  points,  and  the 
crest  tip  passes  the  wall  boundary  without  any  anomaly  developing  in  the  motion  of  points 
near  the  wall  contact  point.  This  run  finally  fails  in  a  way  similar  to  run  E4,  in  which  the 
points  in  the  center  of  the  curl  eventually  become  so  closely  spaced  that  the  simulation 
breaks  down.  The  time  of  failure  is  t  =  50.978  (4152). 

In  run  G3,  which  uses  the  5-point  formula,  the  results  are  nearly  identical  to  run  G2. 
The  time  of  failure  is  f  =  50.977  (4093). 

Another  run  is  done  to  see  if  decreasing  the  number  of  points  added  initially  might 
prevent  the  type  of  failure  in  runs  G2  and  G3  due  to  clustering  of  points  in  the  curl. 

Run  G4  is  identical  in  input  to  run  G3  except  that  point  addition  is  only  applied  once 
at  step  3105,  adding  five  points.  Examining  the  profiles  for  steps  3139  to  3839,  the 
simulation  appears  comparable  to  run  G3  as  the  crest  tip  nears  the  wall  boundary. 
However,  soon  after  the  tip  passes  the  distance  corresponding  to  the  end  of  the  tank,  the 
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lack  of  resolution  in  the  curl  produces  a  failure  most  like  runs  C5  and  C6,  which  had  better 
resolution  in  the  curl  than  run  C4  due  to  the  longer  tank  lengths  used.  This  failure  is  clearly 
shown  in  steps  3849  to  3909.  The  failure  time  is  t  =  50.975  (3939). 

To  test  the  combination  of  the  new  finite  difference  approach  with  another  variation 
of  smoothing,  the  following  run  is  done  using  the  5-point  difference  formula  and  the  same 
smoothing  operation  as  in  run  E2,  in  which  the  wall  point  is  not  used  at  all.  Recall  that  in 
run  E2  the  wall  point  motion  appears  essentially  independent  of  the  adjacent  point  motion, 
and  does  exhibit  an  increase  in  velocity  upon  point  addition. 

In  run  G5,  run  C4  is  again  restarted  at  step  3100  with  smoothing  every  five  steps. 
Point  addition  is  applied  twice,  adding  seven  points  at  3105  and  five  points  at  3110. 
Examining  the  profiles  for  steps  3150  to  3250,  the  jump  in  wall  point  velocity  is  not 
evident;  however,  the  points  closest  to  the  wall  point  begin  to  move  to  the  right  and  the 
point  next  to  the  wall  has  passed  the  wall  point  by  step  3250.  This  abnormal  motion 
continues  until  failure.  The  change  in  computation  of  wall  point  velocity  causes  a  drastic 
difference  in  results  compared  to  run  E2.  The  other  3-  and  4-point  formulas  were  not  tried 
in  combination  with  this  type  of  smoothing,  since  the  5-point  formula  had  worked  well  in 
run  G3. 

Run  G6  attempts  to  solve  the  clustering  problem  observed  in  run  G3  by  using  a 
routine  which  locally  regrids  the  points  in  the  curl  region  when  the  point  spacing  becomes 
less  than  the  smallest  panel  at  the  crest  tip.  The  resulting  run  is  identical  to  run  G3  until 
step  3870,  when  the  local  regridding  routine  is  called.  As  seen  in  the  profiles,  the  points  in 
the  curl  are  evenly  redistributed  using  quadratic  interpolation.  The  run  then  proceeds  until 
an  anomaly  develops  in  the  flow  at  the  points  nearest  the  wall,  similar  to  run  Gl.  The  run 
fails  at  t  -  50.977  (4003).  The  reason  for  this  anomaly  is  not  known,  but  it  should  be  noted 
that  the  failure  occurs  after  the  crest  has  passed  the  position  corresponding  to  the  end  of  the 
tank. 
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The  above  observations  are  summarized  below. 

1.  Runs  G1  through  G3  represent  an  improvement  over  previous  attempts  in  that 
the  increase  in  wall  point  velocity  upon  point  addition  is  no  longer  present. 
Also,  there  is  no  significant  motion  of  wall  points  to  the  left  in  the  later  stage 
of  overturning  as  in  previous  runs.  Although  run  G1  fails  due  to  abnormal 
motion  of  the  points  near  the  wall  after  the  crest  tip  has  passed  the  wall 
boundary,  runs  G2  and  G3  exhibit  no  such  abnormality  and  fail  only  when  the 
points  in  the  center  of  the  curl  become  too  closely  spaced. 

2.  The  use  of  the  new  finite  difference  approach  for  computing  wall  point 
velocity,  in  combination  with  smoothing  which  does  not  directly  alter  the 
present  values  of  the  wall  point  position  and  velocity  potential  but  does  use 
these  values  to  smooth  adjacent  points  appears  to  achieve  a  realistic 
simulation  until  the  crest  tip  has  passed  the  end  of  the  tank.  The  4-  and  5- 
point  formulas  yield  nearly  identical  results. 

3.  Run  G4  shows  again  that  a  minimum  resolution  in  the  curl  region  is  required 
or  the  later  stage  of  overturning  will  not  be  simulated  correctly. 

4.  Run  G5  suggests  that  omission  of  the  wall  point  in  smoothing  is  not  desirable 
when  the  wall  point  velocity  is  computed  using  the  new  finite  difference 
technique.  The  jump  in  wall  point  velocity  is  still  eliminated  as  in  runs  G1  - 
G3,  but  motion  of  the  adjacent  points  becomes  distorted.  This  suggests  that 
coupling  to  the  wall  point  through  smoothing  may  be  a  necessary  constraint 
for  correct  velocity  computation  at  the  adjacent  points. 
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5.3.2  Comparison  of  Successful  Run  with  Experiment 

In  run  C4  in  Section  5.2.2,  it  appeared  that  the  height  of  the  pocket  at  impact  would 
be  about  0.01  providing  the  simulation  could  continue.  From  the  profile  at  step  3793  (i(t)  = 
50.969)  in  run  G3,  the  estimated  pocket  height  at  impact  will  be  closer  to  0.015.  Note  that 
this  height  is  smaller  than  the  initial  length  of  the  free  surface  panels,  0.039.  The  maximum 
crest  height  at  impact  will  be  about  0.19  and  the  wall  run-up  will  be  about  0.16.  Referring 
to  the  comparison  of  run  B6E  with  experiment  in  Section  5.1.2,  the  surface  profile  in  the 
actual  experiment  for  L  =  1 1.667  showed  a  crest  height  of  about  0.21  and  a  wall  run-up  of 
about  0.16  just  prior  to  impact,  and  the  pocket  height  at  impact  is  estimated  from  the  sketch 
to  be  about  0.02  -  0.04.  The  profile  of  the  simulated  wave  as  it  approaches  the  wall  is  thus 
comparable  to  the  actual  experiment,  except  that  the  tank  length  used  in  the  simulation  is  L 
=  1 1 .60  instead  of  1 1 .667. 

Run  G3  is  considered  a  successful  run,  since  curl  resolution  appears  sufficient  and 
there  are  no  apparent  abnormalities  by  the  time  the  crest  tip  reaches  x  =  11.60, 
corresponding  to  the  end  of  the  tank.  This  run  is  therefore  chosen  to  continue  the 
simulation  of  impact  in  the  next  chapter. 

5.4  Summary 

This  effort  to  simulate  wave  breaking  in  a  shorter  tank  has  shown  that  certain 
difficulties  arise  when  overturning  occurs  near  the  vertical  wall  at  the  end  of  the  tank.  In 
the  long  tank,  the  use  of  the  original  smoothing  routine  every  five  time  steps  is  sufficient  to 
allow  the  free  surface  points  to  concentrate  in  the  wave  cusp  and  curl  to  give  good 
resolution  of  the  overturning  process.  As  the  wave  approaches  the  wall  in  the  shorter  tank 
and  begins  to  overturn,  the  resolution  in  the  curl  and  at  the  crest  tip  is  not  as  fine  and  the 
original  points  nearest  the  wall  tend  to  move  to  the  left  into  the  curl  of  the  wave,  causing  a 
further  loss  of  resolution  near  the  wall. 
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In  order  to  adequately  simulate  the  overturning  process,  which  is  characterized  by  an 
upward  velocity  of  the  free  surface  at  the  wall  comparable  to  the  approaching  crest 
horizontal  velocity,  points  are  added  prior  to  overturning  to  increase  resolution  and  both  the 
smoothing  routine  and  velocity  computation  for  the  wall  contact  point  are  modified  to 
eliminate  anomalous  motion  of  the  wall  contact  point  and  adjacent  points.  With  these 
modifications  to  the  program,  a  simulation  with  L  =  1 1 .60  is  finally  obtained  in  which  the 
crest  tip  reaches  the  wall  boundary  position  prior  to  re-entering  the  free  surface  and  without 
any  abnormalities  in  the  flow  up  to  that  point. 

The  problem  of  insufficient  resolution  of  the  free  surface  in  the  cases  of  wave 
breaking  near  the  wall  is  solved  by  the  addition  of  points  prior  to  overturning  and  the 
elimination  of  abnormal  motion  of  the  wall  contact  point  and  adjacent  points.  Note, 
however,  that  these  added  points  eventually  lead  to  failure  because  they  tend  to  cluster 
together  in  the  curl.  In  Run  G3,  this  failure  occurs  well  after  the  crest  tip  has  passed  the  end 
of  the  wall  where  impact  will  be  defined  to  occur,  so  that  simulation  of  the  impact  process 
in  the  next  chapter  should  not  be  affected  by  this  problem.  Run  G6  is  an  attempt  to 
alleviate  the  clustering  problem  by  automatically  regridding  the  points  in  the  curl  region 
when  the  point  spacing  becomes  too  close.  The  initial  regridding  appears  satisfactory,  but 
an  anomaly  eventually  occurs  near  the  wall  which  leads  to  failure. 

The  fact  that  the  addition  of  points  as  described  above  must  be  done  in  a  trial  and 
error  fashion  and  still  can  have  an  undesirable  "side  effect"  of  excessive  clustering  suggests 
the  need  for  a  more  robust  point  distribution  scheme,  especially  in  the  overturning  regime. 
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Chapter  6 

Simulation  of  a  Plunging  Wave  Near  a  Wall:  Impact  Process 


The  idealized  impact  problem  is  posed  numerically  by  creating  a  wall  contact  region 
at  the  crest  at  some  arbitrary  instant  of  time  after  the  tip  of  the  crest  has  passed  the  position 
x  =  L,  as  represented  in  Figure  2-2.  The  essential  changes  to  the  program  in  order  to  model 
impact  are  the  modification  of  existing  subroutines  to  include  the  two  free  surfaces  and  two 
right  wall  contact  regions  after  impact  occurs,  along  with  the  development  of  new 
subroutines  to  define  the  initial  values  of  vertical  position  and  velocity  potential  at  the  new 
free  surface-wall  intersection  points  created  at  impact  of  the  crest  with  the  wall.  Note  that 
in  this  chapter,  the  pressure  in  the  pocket  formed  at  impact  between  the  wave  front  and  the 
wall  is  assumed  to  be  zero;  i.e.,  there  is  no  air  "trapped"  in  the  pocket. 


It  is  important  to  note  that  the  new  wall  contact  points  on  the  free  surfaces  of  the  crest 
are  treated  numerically  in  the  same  way  as  the  original  two  intersection  points  are  prior  to 
impact:  the  values  of  both  velocity  potential  <j)  and  stream  function  i|i  are  specified  in 
computing  influence  coefficients  and  setting  up  the  system  of  equations.  Also,  they  are 
treated  as  Lagrangian  points,  meaning  that  their  vertical  position  and  velocity  potential  are 
obtained  by  integration  of  equations  2.4. 


The  time  at  which  impact  is  defined  to  occur  depends  on  how  milch  of  the  crest  tip  is 
allowed  to  pass  the  postion  x  =  L  before  imposing  the  impact  conditions.  This  can  either  be 
done,  for  instance,  by  specifying  a  minimum  value  of  horizontal  crest  tip  overshoot  after 
which  impact  will  be  imposed,  or  by  specifying  a  maximum  allowable  overshoot  before 
which  impact  must  be  imposed,  or  both.  In  the  present  programs,  a  maximum  allowable 
overshoot  is  specified  and  impact  is  imposed  at  the  first  time  step  at  which  the  crest  tip 
passes  x  =  L  but  is  still  less  than  this  tolerance.  If  the  crest  tip  both  passes  x  =  L  and 
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exceeds  the  tolerance  in  the  same  time  step,  the  simulation  is  restarted  with  time  step  size 
cut  in  half  until  condition  for  impact  is  met. 

Once  the  time  step  at  impact  is  determined,  the  vertical  postion  of  the  crest 
intersection  points  is  set  by  using  quadratic  interpolation  using  values  at  the  point  which  has 
just  passed  x  =  L  as  the  right  endpoint  and  the  two  points  to  its  left  as  the  other  two  nodes. 
The  Jt-position  of  the  intersection  points  is  of  course  x  =  L,  and  the  value  of  y  =  0  is 
assigned  to  the  points.  This  results  in  the  portion  of  the  crest  tip  greater  than  x  -  L  being 
"cut  off.  The  lost  fluid  area  or  volume  is  small  relative  to  the  original  volume,  typically  no 
more  than  0.0001%. 

In  order  to  eventually  compute  pressures  at  the  wall  due  to  crest  impact,  interior 
points  on  the  wall  must  be  added  in  the  crest  contact  region.  This  is  because  the  pressures 
at  the  intersection  points  alone  will  always  be  close  to  zero  if  computed  correctly.  In  one 
version  of  the  impact  program  described  below,  points  are  initially  added  with  even  spacing 
at  the  impact  step,  and  the  spacing  is  set  by  dividing  the  height  of  the  contact  region 
between  the  intersection  points  by  the  average  of  the  lengths  of  the  two  free  surface  panels 
whose  right  endpoints  are  the  intersection  points.  Like  the  points  on  the  wavemaker 
boundary  and  right  wall  boundary,  these  new  points  are  always  evenly  spaced  by  regridding 
the  crest  contact  region  at  every  solution  step  (this  regridding  is  not  be  confused  with  the 
regridding  of  free  surface  points).  The  value  of  \j/  =  0  is  assigned  to  these  points  and  the 
velocity  potential  is  computed  at  each  solution  step. 

Two  methods  of  assigning  initial  values  of  the  velocity  potential  to  the  crest 
intersection  points  are  investigated.  In  the  first  method,  the  velocity  potential  value  is 
obtained  by  quadratic  interpolation  just  as  the  y-position.  In  this  case  the  velocity  potential 
is  assigned  without  any  impact  condition  constraint,  and  the  fact  that  impact  is  occurring  is 
reflected  initially  only  by  the  =  0  values  assigned  to  the  intersection  points  and  to  any 
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added  points  in  the  contact  region.  In  the  second  method,  the  velocity  potential  at  the 
intersection  point  is  computed  subject  to  the  additional  constraint  that  dp/9z  =  0  -  iv  at  that 
point.  This  constraint  is  imposed  by  solving  the  one-sided  3-point  difference  formula 
expressing  5(3/9z  =  u  -  iv  for  <()  at  the  intersection  point  given  the  position  and  complex 
potential  at  the  two  points  closest  to  the  wall  on  the  crest,  the  position  of  the  intersection 
point,  and  the  condition  that  y  =  0  at  the  intersection  point.  This  technique  in  a  sense  is  a 
stronger  imposition  of  the  impact  condition  in  its  redundant  specification  of  zero  normal 
velocity  through  the  initial  value  of  the  velocity  potential  as  well  as  the  condition  y  =  0. 

6.1  Results 

The  results  using  three  versions  of  the  impact  programs  are  presented.  The  primary 
characteristics  of  each  program  are  given  below. 

1 .  Program  1  defines  the  impact  process  only  as  the  creation  of  two  new 
intersection  points.  No  points  are  added  on  the  wall  in  the  impact  region. 

The  purpose  of  this  version  is  to  obtain  an  initial  idea  of  how  the  impact 
simulation  will  work  with  no  other  changes. 

2.  Program  2  is  the  same  as  program  1  but  adds  a  routine  to  check  of  the  free 
surfaces  of  the  crest  at  each  time  step  to  determine  if  any  points  nearest  the 
two  new  intersection  points  have  passed  the  position  corresponding  to  the  end 
of  the  tank.  If  such  a  point  is  found,  a  new  intersection  point  is  defined  at  that 
point  using  one  of  the  two  methods  described  in  the  introduction  and  the  old 
intersection  point  is  eliminated  along  with  any  intervening  points.  This  leads 
to  two  possible  cases,  using  the  top  of  the  crest  as  an  example: 


a.  The  new  intersection  point  is  above  the  old  intersection  point. 
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Physically,  this  corresponds  to  the  initial  contact  region  being 
"enfolded"  by  fluid  which  moves  to  contact  the  boundary  before  the 
initial  contact  region  can  "spread  out".  For  a  spherical  drop  of  fluid 
impacting  a  plane  surface,  Oguz  and  Prosperetti  [8]  give  a  simplified 
analysis  using  mass  conservation  and  the  momentum  equation  to  show 
that  the  initial  line  of  contact  cannot  move  radially  outward  fast 
enough  to  prevent  further  contact  of  the  wall  by  the  fluid  surface, 
unless  the  impact  velocity  is  very  small.  For  a  wave  crest  impacting  a 
wall,  similar  loss  or  destruction  of  the  free  surface  may  well  occur, 
although  the  process  is  clearly  complex. 

b.  The  new  intersection  point  is  below  the  old  intersection  point.  In  this 
case,  this  amounts  to  ignoring  effect  of  fluid  between  the  old  and  new 
intersection  points.  If  this  neglected  fluid  is  part  of  a  thin  jet  which 
has  formed  as  the  initial  impact  region  spreads  out  (the  opposite  case 
of  the  region  being  enfolded),  then  physically  this  may  be  acceptable 
because  the  contribution  to  impact  pressure  in  a  thin  jet  should  be 
small. 

3.  Program  3  is  the  same  as  program  1  but  initially  adds  points  to  the  contact 
region  in  the  impact  step  as  described  in  the  introduction.  The  purpose  of  this 
version  is  to  allow  the  computation  of  impact  pressures  on  the  wall.  For 
simplicity,  the  creation  of  new  intersection  points  after  initial  impact  as  done 
in  program  2  is  not  considered  in  order  to  keep  the  number  of  points  on  the 
wall  constant  during  the  impact  process. 

Each  of  these  programs  is  also  run  both  with  and  without  smoothing  of  free  surface 
points.  The  smoothing  formulas  used  are  identical  to  those  in  the  modified  routine  used  in 
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run  G3  prior  to  impact.  Smoothing  is  applied  over  both  free  surfaces  every  five  time  steps 
after  initial  impact  occurs,  unless  a  subsequent  impact  occurs  (as  in  program  2)  in  which 
case  the  time  step  count  starts  over  and  smoothing  is  delayed  until  five  uninterrupted  steps 
have  passed. 

6.1.1  Summary  of  Attempted  Runs 

In  run  G3  of  Chapter  5,  the  crest  tip  fust  passes  the  wall  at  step  3836  at  t  =  50.972. 
For  the  following  runs,  impact  is  imposed  at  two  different  time  steps  by  restarting  run  G3  at 
steps  3850  (t  =  50.9729)  and  3870  (r  =  50.9738),  respectively.  The  overshoot  tolerance  is 
set  large  enough  so  that  no  backing  up  and  restart  of  the  simulation  is  necessary  and  the  size 
of  the  time  step  is  fixed  at  the  value  before  impact  occurs.  The  impact  time  steps  are  then 
3851  and  3871.  At  3851,  impact  is  imposed  when  the  crest  tip  position  exceeds  L  =  1 1.60 
by  0.0014,  creating  an  initial  contact  length  of  0.0037;  A/  is  set  to  0.000048.  At  step  3871, 
impact  is  imposed  when  the  crest  tip  exceeds  L  =  11.60  by  0.0028,  creating  an  initial 
contact  length  of  0.0054;  At  =  0.000038. 

The  runs  corresponding  to  impact  at  3851  are  designated  by  "H",  and  those 
corresponding  to  impact  at  3871  by  "I".  The  ran  numbers  have  the  following  meaning: 

1  -  Program  1 ,  first  method  of  computing  4>  at  new  points; 

2  -  Program  1,  second  method  of  computing  <|>  at  new  points; 

3  -  Program  2,  fust  method  of  computing  <|>  at  new  points; 

4  -  Program  2,  second  method  of  computing  <j>  at  new  points; 

5  -  Program  3,  first  method  of  computing  <j>  at  new  points; 

6  -  Program  3,  second  method  of  computing  (j)  at  new  points. 

The  further  breakdown  of  run  numbers  refers  to  whether  smoothing  is  used.  For 
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example,  1 1.1  does  not  use  smoothing  and  11.2  does  use  smoothing.  The  run  numbers 
correspond  to  the  free  surface  profiles  given  in  Appendices  H  and  I.  The  abbreviation  "IP" 
will  be  used  below  in  referring  to  the  intersection  points  between  the  wall  and  free  surface. 

In  run  Hl.l,  the  point  next  to  the  top  IP  immediately  overtakes  the  top  IP  at  step 
3852,  leading  to  failure  at  step  3856.  The  initial  velocities  of  both  the  top  and  bottom  IP’s 
are  upward,  29.1  and  19.6,  respectively.  There  is  an  initial  jump  in  the  upward  velocity  at 
the  IP  at  the  bottom  of  the  curl  from  16.5  to  17.9  at  step  3851  when  impact  is  defined,  a 
change  of  $%. 

In  run  H2.1,  the  top  free  surface  points  do  not  tend  to  overtake  the  top  IP  as  in  Hl.l, 
but  the  point  next  to  the  top  IP  moves  rapidly  downward,  deforming  the  crest  and  causing 
failure  at  step  3858.  The  initial  velocities  of  the  top  and  bottom  IP’s  are  49.0  and  24.3, 
respectively.  The  jump  in  velocity  at  the  IP  at  the  bottom  of  the  curl  is  from  16.5  to  18.4,  a 
change  of  12%. 

In  run  H3.1,  the  point  next  to  the  top  IP  which  overtakes  the  IP  in  run  Hl.l  at  step 
3852  is  instead  defined  as  the  new  top  IP  at  the  same  step.  At  step  3851,  the  initial 
velocities  of  the  top  and  bottom  IP’s  are  29.1  and  19.6.  When  the  new  top  IP  is  defined  in 
step  3852,  the  velocity  of  the  top  IP  is  increased  to  31.0  and  the  bottom  IP  assumes  a 
downward  velocity  of  -98.8,  starting  a  downward  jet.  The  top  IP  also  begins  an  upward  jet. 
New  intersection  points  continue  to  be  defined  as  the  simulation  progresses.  Note  that  at 
step  3854  the  point  next  to  the  wall  at  the  bottom  of  the  curl  has  been  disturbed  by  the 
influence  of  the  downward  moving  jet.  Failure  occurs  at  step  3857.  The  initial  velocity 
jump  at  the  IP  at  the  bottom  of  the  curl  is  the  same  as  in  run  Hl.l. 

In  run  H4.1,  the  same  failure  occurs  as  in  run  H2.1  as  the  point  next  to  the  top  IP 
moves  rapidly  downward. 

In  run  H5.1  two  interior  points  are  added  in  the  contact  region  The  profiles  show 
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that  the  free  surface  points  at  the  top  of  the  crest  overtake  the  top  IP,  while  the  free  surface 
points  near  the  wall  at  the  bottom  of  the  crest  form  a  downward  moving  jet.  The  velocity  of 
the  top  IP  is  initially  29.7  and  that  of  the  bottom  point  is  -25.0.  The  simulation  only 
proceeds  for  2  time  steps  before  failure.  At  the  first  impact  step,  3851,  there  is  a  jump  in 
the  upward  velocity  at  the  IP  at  the  bottom  of  the  curl  from  16.5  before  impact  to  17.7  after 
impact  is  imposed,  a  change  of  7%. 

In  run  H6.1,  the  top  free  surface  points  do  not  tend  to  move  past  the  top  IP  and  the 
bottom  free  surface  points  again  form  a  downward  jet.  The  initial  velocities  of  the  top  and 
bottom  IP’s  are  58.4  and  -24.6,  respectively.  This  run  fails  when  the  point  next  to  the  top 
IP  moves  rapidly  downward,  deforming  the  crest,  as  in  runs  H2.1  and  H4.1.  The  initial 
jump  in  velocity  at  the  IP  at  the  bottom  of  the  curl  is  11%. 

The  runs  H1.2,  H2.2,  H3.2,  H4.2,  H5.2  and  H6.2  which  use  smoothing  are  essentially 
the  same  as  their  counterparts  above  which  do  not  use  smoothing  because  of  the  short 
duration  of  the  simulations.  For  this  reason  profile  plots  are  not  given  for  these  runs  in 
Appendix  H. 

In  run  II.  1,  the  point  next  to  the  top  IP  overtakes  the  top  IP  at  step  3875,  causing 
failure.  The  initial  velocities  of  both  the  top  and  bottom  crest  IP’s  are  upward,  1 1 .0  and 
10.5,  respectively,  after  which  they  both  decelerate.  The  jump  in  velocity  at  the  IP  at  the 
bottom  of  the  curl  is  10%. 

Run  11.2  fails  just  as  11.1;  smoothing  does  not  occur  prior  to  failure.  The  initial 
velocities  are  the  same  as  in  1 1 . 1 .  No  plots  are  given  for  this  run. 

Run  12.1  differs  significantly  from  1 1.1  in  that  the  point  next  to  the  top  IP  does  not 
overtake  the  IP,  but  it  does  move  rapidly  downward,  leading  to  failure.  Both  top  and 
bottom  IP’s  initially  move  upward  at  velocities  of  12.8  and  14.2,  and  then  decelerate.  Note 
also  that  the  curl  begins  to  deform  as  clustering  becomes  more  pronounced.  This 
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simultaneous  upward  motion  of  both  IP’s  does  not  seem  physical.  The  jump  in  velcity  at 
bottom  of  the  curl  is  12%. 

Run  12.2  is  very  similar  to  12.1,  but  the  effect  of  smoothing  can  be  seen  in  the  fact 
that  the  point  next  to  the  top  IP  does  not  tend  to  move  downward  as  in  12.1.  Because  no 
restrictions  are  imposed  on  points  other  than  the  initial  IP’s,  the  point  below  the  bottom 
crest  IP  has  passed  x  =  1 1 .6  in  step  3911.  The  initial  velocities  at  impact  are  the  same  as  in 
12.1. 

Run  13.1  is  identical  to  II.  1  up  to  step  3875,  at  which  the  point  which  passes  over  the 
top  IP  is  defined  as  the  new  top  IP.  The  result  is  that  failure  does  not  occur  as  in  1 1.1; 
instead,  the  simulation  proceeds  with  the  bottom  crest  IP  moving  downward  to  form  a  jet 
which  eventually  re-enters  the  free  surface  at  step  3880.  The  top  IP  eventually  forms  an 
upward  jet  as  well.  The  initial  velocities  of  the  top  and  bottom  crest  IP’s  are  1 1.0  and  10.5, 
after  which  the  top  IP  accelerates  and  the  bottom  IP  decelerates.  At  step  3875,  when  the 
new  top  IP  is  defined,  new  top  IP  velocity  jumps  from  21 .6  to  9.8  followed  by  rapid  upward 
acceleration  and  the  bottom  IP  velocity  jumps  from  9.6  to  -3.5,  which  starts  the  downward 
jet.  The  jump  in  velocity  at  the  bottom  of  the  curl  at  impact  is  10%. 

Run  13.2  is  identical  to  run  13.1  until  step  3880,  the  first  step  at  which  smoothing 
occurs;  however,  this  is  also  the  step  at  which  the  bottom  of  the  jet  re-enters  the  free 
surface,  rendering  the  simulation  invalid  there. 

Run  14.1  is  identical  to  ran  12.1  until  step  3893,  when  the  point  next  to  the  top  IP 
which  has  moved  rapidly  downward  eventually  passes  x  =  1 1 .60,  at  which  time  it  is  defined 
as  the  new  top  IP;  however,  it  actually  is  below  the  bottom  IP  and  the  simulation  fails. 

Run  14.2  is  identical  to  run  12.2  until  step  3900,  when  the  point  below  the  bottom 
crest  IP  pass  x  =  1 1 .60  and  is  defined  as  the  new  bottom  IP.  When  this  occurs,  the  top  IP 
velocity  jumps  from  12.9  to  10.7  and  the  velocity  of  the  bottom  IP  jumps  from  8.8  to  3.5, 
after  which  they  both  decelerate. 
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Run  15.1  differs  from  run  1 1.1  in  that  four  interior  points  are  added  in  the  contact 
region  at  impact.  The  results  are  similar  in  that  the  points  next  to  the  top  IP  overtakes  the 
top  IP,  in  this  case  at  step  3876.  An  important  difference,  however,  is  that  the  bottom  IP 
tends  to  form  a  downward  jet,  unlike  run  1 1.1.  The  initial  velocities  of  the  top  and  bottom 
IP’s  at  impact  are  10.6  and  -16.8,  and  the  initial  jump  in  velocity  at  the  IP  at  the  bottom  of 
the  curl  is  10%. 

Run  15.2  is  very  similar  to  run  15.1.  The  application  of  smoothing  at  step  3876  is 
seen  to  delay  the  overtaking  of  the  top  IP  until  one  time  step  later  at  3857.  The  initial 
velocities  of  the  IP’s  are  the  same,  and  the  downward  jet  tends  to  form  in  the  same  way. 

Run  16.1  differs  from  run  12.1  only  in  that  four  interior  points  are  added  in  the  contact 
region  at  impact;  both  runs  use  the  second  method  of  assigning  the  velocity  potential  at  the 
crest  IP’s.  The  difference  in  the  result,  however,  is  dramatic.  This  ran  continues  until  the 
jet  formed  at  the  bottom  of  the  crest  re-enters  the  free  surface  at  the  bottom  of  the  curl.  The 
points  at  the  top  of  the  crest  do  not  form  a  jet  but  do  move  upward  while  decelerating.  The 
initial  velocities  of  the  top  and  bottom  crest  IP’s  are  13.2  and  -14.5,  respectively.  The 
initial  jump  in  velocity  at  the  IP  at  the  bottom  of  the  curl  is  12%.  Compared  to  run  H2,  the 
initial  contact  region  is  larger  but  the  number  of  added  points  is  smaller.  Qualitatively,  this 
simulation  appears  physically  plausible,  in  contrast  to  runs  12  and  14  which  use  the  same 
method  of  assigning  velocity  potential  to  the  IP’s  but  add  no  points  to  the  contact  region  on 
impact. 

Run  16.2  is  identical  to  run  16.1  through  step  3875.  Smoothing  is  applied  at  step 
3876,  but  there  is  little  effect  on  the  result  compared  to  run  16.1.  The  initial  velocities  of 
the  IP’s  are  the  same  as  run  16. 1 . 
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6.1.2  Computation  of  Pressures 

For  run  16.1,  the  pressures  are  computed  at  the  crest  impact  region.  At  the  first  step 
after  impact,  the  initial  pressures  at  the  interior  points  are  negative  with  a  maximum 
magnitude  of  39.9.  The  pressures  at  the  top  and  bottom  intersection  points  at  this  step  are 
-0.15  and  -6.6.  The  pressures  at  the  interior  points  remain  negative  but  decrease  until  step 
3874,  where  there  is  a  positive  value  of  29.9.  For  all  subsequent  steps  the  interior  pressures 
are  positive,  attaining  a  maximum  value  of  48.9  at  step  3876. 

In  order  to  obtain  more  imformation  about  the  computed  impact  pressure,  a  series  of 
variations  on  run  16.1  is  done  in  which  the  time  step  size  and  number  of  points  in  the 
contact  region  are  changed.  The  time  step  sizes  used  are  0.00004,  0.00001  and  0.000001. 
At  each  time  step  size,  the  number  of  points  added  are  either  two  or  five.  Plots  showing  the 
time  history  of  the  pressure  profile  on  the  impact  region  of  the  wall  for  the  case  of  At  = 
0.00001  are  shown  in  Appendix  J.  It  is  found  that  the  change  in  time  step  size  has 
negligible  effect  on  the  results  for  a  given  number  of  points;  however,  changing  the  number 
of  points  has  a  marked  difference  on  the  pressures  computed  at  a  given  time  step.  For  the 
case  of  four  total  points  (two  added)  in  the  contact  region,  although  the  pressure  values  are 
initially  negative  they  eventually  turn  positive  with  a  maximum  of  about  50.  When  five 
points  are  added,  the  initial  negative  pressures  are  almost  three  times  as  large.  The  pressure 
values  never  become  completely  positive  as  they  do  in  the  case  with  fewer  points,  and  the 
maximum  positive  value  attained  is  about  one  third  as  large. 

It  should  be  noted  that  the  pressure  computations  are  performed  after  the  simulation 
using  a  separate  analysis  program;  there  is  a  possibility  that  the  simulations  could  be  valid 
and  the  analysis  inaccurate  However,  the  same  basic  computation  of  pressure  is  done  for 
the  impact  region  as  is  used  in  computing  the  work  input  from  the  wavemaker,  which  has 
been  reasonably  well-tested.  Examining  the  components  of  the  pressure  calculation  show 
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that  the  large  negative  values  are  associated  with  the  term  Dty/Dt  in  Bernoulli’s  equation 
which  is  computed  using  a  three-  point  difference  in  time. 

6.1.3  Comparison  with  Experiment 

In  the  experiment  of  Chan,  the  typical  range  of  peak  impact  pressures  measured  is 
given  as  5pC 2  -  10pC2  where  p  is  the  water  density  and  C  is  the  linear  wave  phase  speed 
corresponding  to  the  center  frequency  of  the  wavemaker  velocity  spectrum,  or  1.70  m/s 
[1,2],  The  largest  impact  pressure  measure  was  21pC2.  In  terms  of  the  non- 
dimensionalization  described  at  the  beginning  of  Chapter  2,  the  measured  non-dimensional 
impact  pressures  are  in  a  range  from  2.5  -  5  with  the  largest  value  being  about  10.  No  large 
negative  pressures  were  noted  in  the  experiment. 

The  positive  pressures  computed  for  the  simulation  are  of  the  same  order  of 
magnitude  as  those  measured  in  the  experiment,  but  there  is  as  yet  no  explanation  for  the 
computed  negative  pressures.  The  duration  of  impact  in  the  simulation  is  about  0.0003, 
defined  from  impact  until  the  jet  at  the  bottom  of  the  crest  re-enters  the  free-surface. 

Recall  that  air  compressibility  in  the  pocket  has  not  yet  been  included  in  the  model. 
Another  unknown  factor  is  the  effect  of  using  a  shorter  tank  length  than  the  experiment  in 
order  to  get  impact  to  occur. 

6.2  Summary 

The  following  observations  are  made  concerning  the  techniques  tried  in  the 
simulation  of  wave  impact  on  a  vertical  wall: 

1.  Although  the  use  of  the  second  method  of  computing  velocity  potential  at  the 
new  crest  intersection  points  tends  to  prevent  adjacent  points  from  overtaking 
them,  the  resulting  flow  when  no  points  are  added  to  the  interior  of  the  crest 
contact  region  at  impact  appears  non-physical. 
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2.  When  the  first  method  of  computing  velocity  potential  at  the  intersection 
points  is  used,  coupled  with  the  definition  of  new  intersection  points  as 
adjacent  points  overtake  the  original  points  (e.g.,  runs  13.1  and  13.2),  the  flow 
appears  plausible  even  though  no  interior  points  are  added  at  impact.  A 
downward  jet  forms  at  the  bottom  of  the  crest  which  eventually  re-enters  the 
free  surface,  while  the  free  surface  at  the  top  of  the  crest  accelerates  upward. 

3.  When  the  second  method  of  computing  velocity  potential  is  used  coupled 
with  the  addition  of  points  in  the  contact  region  at  impact,  the  flow  again 
appears  plausible  (run  16.1).  If  the  first  method  is  used,  then  the  adjacent 
points  will  overtake  the  new  intersection  points.  This  could  be  remedied  by 
defining  these  overtaking  points  to  become  new  intersection  points  when  they 
reach  the  wall  boundary  position,  as  use  of  program  2  did  in  runs  13.1  and 
13.2.  This  was  not  attempted  due  to  the  desire  to  keep  the  number  of  points 
on  the  wall  constant  for  the  pressure  computations. 

Note  that  the  "H"  runs  which  use  the  second  method  of  computing  velocity  potential 
at  new  intersection  points  fail  due  to  the  rapid  downward  motion  of  the  point  next  to  the  top 
crest  intersection  point.  This  could  probably  be  improved  if  smoothing  were  applied  more 
often,  such  as  every  time  step,  since  similar  motion  of  this  same  point  is  "corrected"  by 
smoothing  in  run  12.2. 

The  computations  of  pressure  show  results  independent  of  time  step  size  but  highly 
dependent  on  the  number  of  points  added  at  impact  to  the  contact  region.  The  initial 
pressures  computed  are  negative  in  all  test  cases,  and  in  general  grow  more  positive  with 
time.  The  highest  positive  pressure  computed  is  about  50  for  the  case  with  only  two  points 
added,  while  the  negative  pressures  of  largest  magnitude  occur  in  the  case  of  five  points 
added.  The  positive  pressures  obtained  are  comparable  in  order  of  magnitude  to  the 
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experimental  values,  but  there  is  no  explanation  as  yet  for  the  initial  negative  pressures. 
These  simulations  do  not  include  the  compression  of  air  in  the  pocket  formed  by  the 
impacting  wave. 
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Chapter  7 

Investigation  of  Criteria  for  a  Robust  Point  Distribution  Scheme 

The  results  of  simulations  of  plunging  breaking  waves  in  Chapter  4  and  Chapter  5 
have  shown  that  the  distribution  of  points  on  the  free  surface  has  a  major  effect  on  the 
outcome  of  a  simulation.  "Distribution"  is  meant  here  in  a  general  sense,  referring  to  the 
variation  of  both  the  position  and  complex  potential  from  point  to  point.  The  distribution  of 
points  is  initially  given  at  t  =  0  by  setting  the  velocity  potential  at  all  points  equal  to  zero 
and  using  a  total  number  of  points  so  that  the  initial  panel  density  is  about  40  panels  per 
wavelength  based  on  a  linear  wave  with  a  frequency  ofto  =  2.0,  the  frequency  below  which 
most  of  the  input  energy  is  concentrated.  As  the  simulation  proceeds,  the  point  distribution 
obviously  changes  in  time.  The  regridding,  smoothing  and  point  addition  routines  used 
thus  far  change  the  point  distribution  in  some  way  when  they  are  applied  at  a  particular  time 
step.  The  regridding  routine  has  been  used  as  a  method  of  suppressing  sawtooth 
instabilities  by  intermittently  redistributing  points  more  or  less  evenly  using  quadratic 
interpolation  with  the  linear  distance  between  points  as  a  parameter.  Consequently, 
regridding  is  not  desired  when  better  resolution  of  the  free  surface  is  required  in  the 
overturning  phase.  Smoothing  is  then  substituted  for  regridding  to  capture  the  overturning 
event  by  allowing  the  concentration  of  Lagrangian  points  while  still  serving  to  suppress 
sawtooth  instabilities.  Smoothing  is  in  a  sense  a  local  redistribution  of  points  which  finds  a 
mean  quadratic  curve  through  five  points  and  then  redefines  the  points  to  be  on  this  curve. 
In  the  runs  using  the  short  tank  lengths,  the  need  for  even  better  resolution  in  the 
overturning  phase  than  that  provided  by  the  existing  points  is  met  by  using  a  routine  which 
redistributes  points  by  adding  points  based  on  quadratic  interpolation  when  point  spacing 
exceeds  an  arbitrary  multiple  of  the  minimum  panel  size  at  the  crest  tip. 
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The  application  of  these  point  distribution  routines  in  the  simulations  discussed  thus 
far  has  been  essentially  a  trial  and  error  process.  In  order  to  develop  a  program  which  can 
ultimately  simulate  wave  breaking  and  impact  without  such  an  ad  hoc  use  of  separate  point 
distribution  techniques,  development  of  a  more  robust  point  distribution  scheme  is 
necessary  which  can  automatically  control  the  distribution  of  points  to  optimize  the 
simulation  from  the  standpoint  of  computational  time  and  accuracy.  Such  a  scheme  would 
need  certain  criteria  for  distributing  points  based  upon  comparison  of  selected  parameters 
with  critical  values  of  these  parameters.  It  could  conceivably  use  the  existing  point 
distribution  routines  in  a  more  sophisticated  way  or  could  even  be  fundamentally  different 
from  them. 

The  objective  of  this  chapter  is  to  examine  possible  parameters  and  types  of  criteria 
which  may  eventually  prove  useful  in  such  a  robust  distribution  scheme.  The  approach 
taken  is  to  examine  representative  cases  where  simulations  have  failed  in  an  attempt  to  find 
quantities  which  exhibit  marked  changes  near  the  failure  point.  Such  quantities  may  then 
be  candidates  for  use  as  parameters  in  defining  the  criteria.  The  types  of  observed  failures 
are  summarized  below. 

1 .  Without  the  application  of  any  smoothing  or  regridding  routine  at  all,  a 
simulation  will  fail  due  to  a  "sawtooth"  type  of  instability  in  the  free  surface. 

In  run  B4  the  simulation  is  started  at  t  =  0  without  regridding  or  smoothing, 
and  sawtooth  instability  develops  near  the  wavemaker  which  causes  failure;  in 
run  D8  smoothing  is  stopped  just  before  overturning  and  a  sawtooth 
instability  forms  in  the  curl  of  the  wave  and  near  the  wall  which  leads  to 
failure. 

2.  When  regridding  alone  is  used  from  the  start  of  the  simulation,  failure  occurs 
as  a  wave  crest  steepens  and  tends  to  overturn.  Such  failures  occur  at  the 


-75- 


peak  of  the  crest  and  seem  to  suggest  insufficient  resolution  of  the  free 
surface  in  this  region  due  to  regridding.  Run  B1  is  an  example. 

3.  In  simulations  where  smoothing  is  being  applied  and  overturning  is 

I  progressing,  failure  can  occur  due  to  insufficient  resolution  in  the  curl  of  the 

wave  and  near  the  wall.  Runs  B6E  and  C4  are  examples. 

4.  In  simulations  where  smoothing  is  being  applied  and  overturning  is 
progressing,  point  addition  in  the  curl  region  to  improve  resolution  can 
eventually  lead  to  failure  due  to  "excessive"  resolution  as  points  at  the  bottom 
of  the  curl  tend  to  be  convected  upward  into  the  center  of  the  curl.  Failure 
occurs  when  these  points  become  too  closely  spaced.  Runs  E2  and  G3  are 
examples. 

Recall  in  both  Chapters  4  and  5  that  several  simulations  "fail"  due  to  an  apparent 
"weak  breaking"  of  the  peak  preceeding  the  one  which  is  known  to  overturn  in  experiment. 
In  physical  terms  these  runs  are  not  necessarily  failures;  they  could  be  accurate  at  least  in 
that  the  first  peak  may  actually  have  experienced  a  degree  of  breaking;  this  cannot  be 
verified  without  repeating  the  experiment.  In  order  to  permit  the  simulation  to  proceed  to 
the  overturning  of  the  second  peak,  a  combination  of  regridding  and  smoothing  is  used  to 
suppress  the  overturning  of  the  first  peak.  The  point  is  that  the  use  of  regridding  plus 
smoothing  in  these  cases  is  an  artificial  means  of  achieving  an  end,  and  in  general  may  be 
of  limited  use  trying  to  simulate  problems  without  any  a  priori  knowledge  of  the  outcome. 
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7.1  Selection  of  Parameters  and  Test  Cases 

All  of  the  observed  failure  types  described  above  are  localized  in  some  way.  This 
fact  suggests  that  any  criteria  which  will  be  used  to  determine  point  distribution  on  the  free 
surface  should  be  based  upon  parameters  related  to  local  flow  properties  and  local  point 
distribution  rather  than  on  global  quantities.  Three  such  parameters  are  defined  below. 

One  parameter  to  be  examined  is  designated  SAW  and  is  intended  to  provide  an 
indication  of  the  presence  of  sawtooth  instability  on  the  free  surface.  It  is  defined  as 
follows: 

N 

SAW=  £(-l)"(Aa) 

n=i 

where  the  free  surface  panels  are  numbered  from  1  to  A  and  Act  is  the  relative  radian  angle 
between  a  panel  and  the  previous  panel  (positive  counterclockwise)  divided  by  tc.  The  idea 
is  that  a  "sawtooth"  geometry  of  the  panels  will  cause  an  alternating  sign  on  Aa  from  panel 
to  panel,  which  will  tend  to  increase  the  magnitude  of  SAW. 


The  second  parameter  to  be  examined  is  a  non-dimensional  estimate  of  local 
curvature,  designated  CUR,  and  is  defined  at  each  panel  vertex  as  follows: 

CUR  =  ~  Sq 


where  S  is  the  sum  of  the  half-panel  lengths  on  either  side  of  the  vertex  and  S0  is  the 
original  panel  length  at  t  =  0. 

The  third  parameter  is  a  measure  of  the  rate  of  contraction  or  expansion  of  a  panel 
relative  to  the  average  velocity  of  the  panel  endpoints.  It  is  designated  CON  and  is  defined 
as  follows: 
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M  1 

CON  = 

At  V 


where  A /  is  the  change  in  panel  length  and  V  is  the  average  velocity  magnitude  of  the 
endpoints  of  the  panel. 

There  are  essentially  two  approaches  to  examining  such  parameters.  One  is  to  track 
the  value  of  the  parameter,  or  perhaps  its  maximum  or  minimum,  in  time.  The  other  is  to 
choose  a  particular  time  step  and  look  how  the  value  of  the  parameter  varies  spatially  along 
the  free  surface,  assuming  it  is  defined  at  each  free  surface  point  or  panel. 

The  runs  selected  as  representative  cases  are  runs  Bl,  B4,  C4  and  G3. 

7.2  Definition  of  Criteria 

Once  parameters  of  interest  have  been  identified,  their  correlation  with  the  failure  of 
a  simulation  must  be  determined.  If  a  parameter  value  changes  markedly  as  failure 
approaches,  then  it  may  be  related  to  the  reason  for  failure.  This  is  especially  true  if  the 
values  of  the  same  parameter  in  a  run  which  does  not  fail  at  the  same  point  are  considerably 
different  from  the  values  associated  with  failure. 

If  a  strong  correlation  does  exist  between  a  parameter  and  failure,  then  critical  values 
of  the  parameter  can  conceivably  be  found  to  define  a  criterion  for  the  redistribution  of 
points  in  such  a  way  as  to  prevent  failure.  The  robust  point  distribution  algorithm  would 
control  the  starting,  stopping,  frequency  and  method  of  redistribution  based  upon  such 
criteria.  Ideally,  it  should  cause  minimal  distortion  of  the  solution  while  preventing 
instabilities. 

In  the  following  examples  using  the  runs  cited  above,  simple  criteria  are  suggested 
based  upon  the  three  parameters  described  which  could  improve  the  use  of  the  existing 
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point  distribution  routines  by  making  their  application  more  automated  and  rational.  The 
actual  implementation  of  these  suggested  criteria  has  not  yet  been  attempted. 

In  the  first  example,  the  failure  due  to  sawtooth  instablity  is  addressed.  The  values  of 
SAW  are  compared  between  run  Bl,  which  shows  no  sawtooth  instabilities,  and  run  B4, 
which  does  develop  such  instability  near  the  wavemaker. 

Figure  7-1  shows  time  histories  of  SAW  for  the  two  runs.  Note  that  around  t~  27  the 
value  of  SAW  in  run  B4  becomes  significantly  larger  than  the  value  in  run  Bl,  exceeding 
0.1.  Examination  of  the  free  surface  profiles  for  run  B4  in  Appendix  B  shows  that  this  is 
also  roughly  the  time  when  the  sawtooth  instability  becomes  apparent  in  the  plots  (i.e.,  at 
step  1000).  This  suggests  that  SAW  could  be  used  in  a  criterion  which  will  start  smoothing 
or  regridding  only  when  the  value  of  SAW  exceeds,  for  instance,  0.1,  and  stop  once  the 
value  falls  below  0.1  again. 

In  the  second  example,  failure  due  to  lack  of  resolution  associated  with  regridding  is 
considered  by  examining  values  of  the  curvature  parameter  CUR  in  run  B 1 .  The  issue  here 
is  whether  a  criterion  can  be  found  which  can  control  when  to  stop  regridding  and  start 
smoothing  so  that  sufficient  resolution  of  an  overturning  event  can  be  achieved.  It  is 
possible  that  regridding  may  fail  when  the  curvature  of  the  steepening  crest  becomes  too 
large.  About  100  steps  prior  to  failure  at  the  first  peak,  the  value  of  CUR  at  the  points 
defining  the  peak  is  about  two  orders  of  magnitude  larger  than  at  neighboring  points.  At 
step  2636,  for  example,  the  variation  in  CUR  at  the  first  peak  is  from  -0.0067  to  -0.455  and 
then  back  to  0.0061  over  a  span  of  eight  points.  This  suggests  that  a  criterion  could  be 
defined  which  locates  a  peak  and  stops  regridding  and  starts  smoothing  when  the  magnitude 
of  CUR  exceeds  a  critical  value  at  the  peak,  such  as  0.01. 

In  the  third  example,  the  type  of  failure  due  to  lack  of  resolution  in  the  curl  of  an 
overturning  wave  near  a  wall  is  addressed.  Run  C4,  which  fails  due  to  such  lack  of 
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resolution,  is  compared  with  run  G3,  which  proceeds  further  than  C4  because  of  the 
addition  of  points  prior  to  overturning.  The  main  difference  here  is  the  value  of  the 
parameter  CON  at  panels  near  the  wall.  In  run  C4,  at  step  3540  ( t  =  50.962)  the  values  of 
CON  at  the  three  panels  closest  to  the  wall  from  left  to  right  are  -0.41,  -3.37  and  2.60.  At 
this  step  overturning  is  in  progress  (see  profiles  in  Appendix  C).  The  value  of  2.60  for  the 
panel  next  to  the  wall  indicates  that  it  is  expanding  rapidly  relative  to  other  panels,  while 
the  negative  values  indicate  contraction  of  the  adjacent  panels.  This  is  consistent  with  the 
observed  motion  to  the  left  of  the  point  next  to  the  wall  in  the  profiles.  In  contrast,  for  run 
G3  at  step  3693  (r  =  50.963)  the  corresponding  values  of  CON  are  -0.06,  -0.07  and  0.03. 
This  suggests  that  a  criterion  might  be  developed  which  would  add  points  to  increase 
resolution  based  on  the  increasing  magnitude  of  CON  at  panels  near  the  wall.  It  is  not  clear 
what  a  critical  value  of  this  parameter  might  be. 

The  fourth  type  of  observed  failure  which  occurs  when  points  in  the  curl  cluster 
together  has  not  yet  been  linked  to  a  parameter.  However,  the  method  of  local  reridding 
used  in  run  G6  to  alleviate  the  clustering  problem  observed  in  run  G3  still  may  be  worth 
pursuing  as  a  means  of  redistribution  once  an  improved  criterion  for  invoking  it  is 
developed. 

7.3  Summary 

Although  a  robust  point  distribution  scheme  has  yet  to  be  developed,  some  non- 
dimesional  parameters  related  to  local  flow  characteritics  have  been  identified  which  are 
correlated  with  various  observed  types  of  failure.  Criteria  based  upon  these  parameters 
have  been  suggested  which  may  improve  the  application  of  the  existing  types  of 
distribution  routines  and  thereby  lay  the  groundwork  for  the  desired  robust  scheme. 
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Chapter  8 
Conclusions 


8.1  Summary  of  Work 

This  research  attempts  to  extend  the  use  of  fully  nonlinear  potential  flow  theory  from 
the  simulation  of  a  deep-water  plunging  breaking  wave  by  Dommermuth  et  al.  [3]  to  the 
simulation  of  a  deep-water  plunging  breaking  wave  impacting  a  vertical  wall.  The 
experimental  basis  for  this  numerical  study  is  the  work  of  Chan  [1,2]  dene  on  the 
measurement  of  impact  pressures  and  kinematics  of  a  deep-water  plunging  wave  impacting 
a  vertical  wall.  The  numerical  wavemaker  velocity  input  used  is  the  same  as  that  used  by 
Dommermuth  et  al.,  only  the  numerical  tank  is  shortened  as  the  actual  tank  was  in  the 
experiment  in  order  to  force  overturning  of  the  wave  to  occur  near  the  wall  so  that  wall 
impact  is  possible. 

8.1.1  Simulations  for  Long  Tank 

The  simulation  of  a  deep  water  plunging  wave  in  a  long  tank  of  L  =  20  is  repeated  in 
order  to  investigate  the  dependence  of  results  upon  the  use  of  smoothing,  regridding  or 
both.  It  is  confirmed  that  sawtooth  instabilities  develop  on  the  free  surface  in  this  case 
when  the  simulation  is  started  without  the  use  of  any  smoothing  or  regridding.  The 
simulation  is  complicated  by  the  apparent  "weak  breaking"  tendency  of  the  peak  preceding 
the  peak  which  eventually  overturns  and  plunges.  This  difficulty  is  overcome  by  the  use  of 
regridding  plus  smoothing  for  a  limited  interval  which  suppresses  the  breaking  of  the  first 
peak,  after  which  smoothing  alone  is  used  which  allows  simulation  of  the  overturning  of  the 
second  peak  to  continue. 
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,  The  results  of  the  simulations  are  found  to  be  most  dependent  upon  the  ending  time 

of  the  regridding  plus  smoothing  interval.  As  the  time  of  starting  smoothing  alone 
decreases,  the  time  at  which  re-entry  occurs  is  observed  to  increase.  For  the  simulations 
which  continue  to  re-entry,  the  maximum  difference  in  time  of  re-entry  is  0.15,  or  about 
0.3%  of  the  average  total  simulation  time.  The  maximum  deviation  in  re-entry  x-position  is 
0.2,  or  about  1.7%  of  the  average  x-position  at  re-entry. 

i  The  relative  global  conservation  of  energy  is  also  dependent  on  the  timing  of 

regridding  and  smoothing  application.  The  runs  with  the  longest  and  shortest  re-entry  times 
are  found  to  have  the  worst  conservation  of  energy,  gaining  about  7%  in  fluid  energy 
relative  to  total  wavemaker  work  over  the  duration  of  the  simulation. 

8.1,2  Simulations  for  Short  Tank:  Pre-impact 

The  same  time  sequence  of  regridding,  regridding  plus  smoothing  and  smoothing 
alone  that  is  used  in  run  A3B  to  successfully  simulate  a  plunging  breaker  in  the  long  tank  is 
applied  to  the  short  tank  simulations  after  noting  that  the  difficulty  with  the  weak  breaking 
tendency  of  the  first  peak  still  exists.  In  runs  where  no  regridding  and  smoothing  is  used, 
the  development  of  a  sawtooth  instability  of  the  free  surface  is  observed. 

As  in  the  simulation  for  the  tank  of  L  =  20,  breaking  in  the  simulations  using  a 
shorter  tank  occurs  sooner  than  in  the  actual  experiment.  Initially,  a  length  of  L  =  1 1.667  is 
used  to  coincide  with  one  of  the  cases  of  wall  impact  investigated  by  Chan  [1,2],  but 
overturning  and  near  re-entry  in  the  successful  simulation  occurs  before  the  crest  tip 
reaches  the  horizontal  postion  corresponding  to  the  end  of  the  tank.  After  experimenting 
with  several  shorter  tank  lengths,  a  simulation  is  achieved  with  L  =  1 1 .60  in  which  the 
overturning  crest  appears  as  though  it  will  reach  the  end  of  the  tank,  but  failure  occurs  due 
to  lack  of  resolution  of  the  free  surface  near  the  wall. 
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In  order  to  finally  obtain  a  simulation  with  L  =  1 1.60  which  continues  to  the  point  of 
impact  without  any  apparent  abnormalities  in  the  flow,  the  following  techniques  are  used: 

1.  Points  are  added  in  the  curl  region  and  near  the  wall  prior  to  the  start  of 
overturning  to  increase  resolution. 

2.  The  smoothing  of  points  near  the  wall  is  changed  so  that  the  position  and 
velocity  potential  at  the  intersection  point  are  not  smoothed  themselves  but 
are  used  in  smoothing  the  adjacent  two  points.  The  original  method  smoothed 
the  wall  point  values  as  well. 

3.  The  computation  of  vertical  velocity  at  the  intersection  point  is  modified  from 
a  one-sided  3-point  difference  of  complex  potential  along  the  free  surface  to  a 
one-sided  5-point  difference  of  velocity  potential  along  the  right  wall. 

Important  observations  from  these  simulations  are  the  following: 

1 .  There  is  a  minimum  degree  of  resolution  required  in  the  curl  region  to 
successfully  simulate  the  overturning  process. 

2.  Point  addition  prior  to  overturning  can  cause  difficulty  later  in  the  overturning 
process  as  points  tend  to  cluster  in  the  center  of  the  curl.  When  spacing 
becomes  too  close  the  simulation  breaks  down. 

3.  The  modified  finite  difference  scheme  for  computing  intersection  point 
velocity  coupled  with  the  modified  smoothing  routine  eliminates  the  problems 
of  the  jump  in  intersection  point  velocity  upon  point  addition  and  anomalous 
motion  of  points  nearest  the  wall.  In  this  case,  the  3-point  finite  difference 
scheme  eventually  allows  abnormal  flow  near  the  wall,  but  the  4-  and  5-point 
schemes  appear  sufficient  and  give  nearly  identical  results. 


At  the  point  where  the  crest  tip  reaches  the  wall,  the  crest  height,  wall  run-up  and  curl 
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dimensions  for  the  simulated  wave  with  L  =  1 1 .60  are  comparable  to  the  experiment  for  the 
case  where  L  =  1 1.667. 

8.1.3  Simulations  for  Short  Tank:  Post  Impact 

The  idealized  impact  problem  is  initially  posed  by  creating  two  new  wall-free  surface 
intersection  points  at  the  crest  after  the  crest  has  passed  the  horizontal  position 
corresponding  to  the  end  of  the  tank.  Two  techniques  are  examined  for  assigning  the  initial 
velocity  potential  to  these  new  intersection  points,  and  a  technique  is  tried  for  defining  new 
intersection  points  as  other  free  surface  points  later  reach  the  wall  boundary.  The  effect  of 
adding  interior  points  to  the  contact  region  at  impact  is  also  investigated,  and  such  point 
addition  is  observed  to  produce  the  most  plausible  flows  in  a  qualitative  sense. 

The  preliminary  results  of  pressure  computations  for  the  impact  region  without  the 
effect  of  air  compressibility  are  puzzling.  The  initial  pressures  are  negative,  and  tend  to 
become  more  positive  with  time.  Also,  the  magnitude  of  the  initial  pressures  is  seen  to  be 
highly  dependent  on  the  number  of  points  added  to  the  contact  region  at  impact.  Although 
positive  pressures  are  eventually  obtained  which  are  of  the  same  order  of  magnitude  as 
those  measured  in  the  experiment,  there  is  as  yet  no  explanation  for  the  computed  negative 
pressures. 

8.2  Recommendations  for  Future  Work 

8.2.1  Point  Distribution  Schemes 

The  further  investigation  of  parameters  which  may  be  correlated  to  observed  types  of 
failure  should  be  continued  with  the  objective  of  improving  the  methods  of  distributing 
points  on  the  free  surface  and  ultimately  developing  a  robust  algorithm  which  can  ensure 
effective  simulations  for  a  variety  of  nonlinear  water  wave  problems.  In  particular,  the 
criteria  suggested  in  Chapter  7  could  be  implemented  to  determine  their  effectiveness. 
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8.2.2  Impact  Problem 

The  idealized  impact  problem  using  potential  flow  is  observed  in  some  cases  to  yield 
qualitatively  plausible  flows.  The  preliminary  computations  of  pressure,  however,  are  not 
easily  accepted  as  physical.  The  negative  pressures  obtained  on  impact  could  be  the  result 
of  the  following: 

1.  An  error  in  the  main  program  which  solves  the  BVP  and  integrates  the  free 
surface  boundary  conditions. 

2.  An  error  in  the  analysis  program  which  computes  the  pressures  from  the  main 
program  output. 

3.  Basic  assumptions  concerning  how  the  initial  impact  condition  is  to  be 
imposed;  for  example,  the  method  of  initially  assigning  velocity  potentials  to 
the  new  intersection  points  and  the  discontinuity  introduced  by  "cutting  off’ 
an  arbitrary  amount  of  the  crest  tip  to  define  the  impact  region  Note  that  a 
value  of  y  =  0  is  assigned  to  the  points  on  the  wall  in  the  impact  region;  other 
constant  values  of  y  should  be  tried  to  see  if  there  is  any  affect  on  the  results. 

In  addition  to  investigating  the  above  issues,  the  following  suggestions  are  made: 

1.  Include  air  compressibility  in  the  simulation  by  modifying  the  dynamic  free 
surface  boundary  condition  in  the  "pocket"  to  impose  a  pressure  based  upon 
the  volume  of  air  in  the  pocket  and  an  ideal  gas  law,  perhaps  adiabatic.  This 
is  an  essential  next  step  if  results  are  to  be  carefully  compared  with  those  of 
Chan  [1,2],  who  surmised  that  measured  pressure  oscillations  during  impact 
were  due  to  the  dynamics  of  the  air  trapped  at  impact. 

2.  Investigate  possible  point  addition  techniques  and  time  step  control 
particularly  suited  to  the  impact  regime,  where  high  velocities  are  computed 
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for  the  intersection  points  created  at  impact.  In  the  computations  done  so  far, 
the  time  step  size  has  been  held  constant  throughout  impact. 


3.  Investigate  the  possible 
experimental  results. 


scaling  of  the  computational  problem  to  match 
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Appendix  A 

Refer  to  Chapter  4  for  discussion. 
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Refer  to  Chapter  5  for  discussion. 
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Refer  to  Chapter  5  for  discussion. 
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Refer  to  Chapter  5  for  discussion. 
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Refer  to  Chapter  5  for  discussion. 
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Refer  to  Chapter  5  for  discussion. 
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Refer  to  Chapter  6  for  discussion. 
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Refer  to  Chapter  6  for  discussion. 
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115-10  115.20  115.30  115.40  115.50  115.61 


A 


50 


A 


115  10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115  90  116.00  116.10 


A 


115.10  115.20  115.30  115  40  115.50  115.60  115.70  115.80  115.90  116.00 


5.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


‘115.10  115.20  115.30  115  40  115.50  115  60  115.70  115. BO  115.90  116  00  116  10 


115.10  115.20  115.30  115.40  115.51 


I 


I 


» 


► 

► 


15.2  3872  50.9739 


K 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


115.10  115.20  115.30  115.40  115.51 


A 


5.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


5.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90 
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*115. 10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  11600  116.10 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115  90  116.00  116.10 


A 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


115.10  115.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116.00  116.10 


53 


5.20  115.30  115.40  115.50  115.60  115.70  115.80  115.90  116  00  116.10 


535 


115.10  115.20  115.30  115.40  115.50  115.60  115  70  115.80  115.90  116.00  116.10 
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Appendix  J 


Refer  to  Chapter  6  for  discussion. 


Time  history  of  impact  pressure;  time  step  size  =  0.00001,  no.  of  points 


Time  history  of  impact  pressure;  time  step  size  =  0.00001,  no.  of  points 


Time  history  of  impact  pressure;  time  step  size  =  0.00001,  no.  of  points 
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